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

Abstract

The application discloses a three-dimensional anisotropic medium space-time domain seismic forward method, device, medium and product, relates to the technical field of forward simulation of geophysical exploration, and comprises the following steps: based on a preset improved time high-order difference operator and a preset improved space high-order difference operator, respectively discretizing time derivatives and space derivatives in a three-dimensional VTI medium acoustic wave equation; based on the three-dimensional VTI medium acoustic wave equation, simplifying time-discrete functions and space-discrete functions to obtain a simplified space-time domain dispersion relation; based on a target optimization algorithm, solving the simplified space-time domain dispersion relation to obtain optimized time item high-order difference coefficients and optimized space item high-order difference coefficients; and based on the optimized time item high-order difference coefficients, the optimized space item high-order difference coefficients, the time-discrete functions and the space-discrete functions, extrapolating and solving the three-dimensional VTI medium acoustic wave equation. The application can effectively improve the numerical precision in anisotropic wave field simulation.

Description

Technical Field

[0001] The present application relates to the field of geophysical exploration forward modeling technology, 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, the exploration and development of established oil fields places higher demands on the accuracy and efficiency of geophysical exploration technology. Traditional seismic exploration technology is primarily based on the assumption of isotropic media, leaving significant room for improvement in accuracy. With the increasing complexity of exploration targets and the continuous advancement of exploration technology, the development of new seismic exploration theories, methods, and technologies based on anisotropy theory, even for more complex media characteristics, is of great significance. Seismic data processing plays a connecting role in seismic exploration. Among them, seismic migration imaging and waveform inversion can effectively help understand subsurface structural or structural characteristics. Imaging and inversion methods and techniques rely on high-precision wavefield continuation strategies. Therefore, it is necessary to develop high-precision seismic wavefield numerical simulation methods and technologies to complement imaging and inversion techniques based on anisotropy theory.

[0003] To date, researchers have developed a variety of numerical algorithms for the numerical solution of wave equations, such as pseudospectral methods, finite element methods, and finite difference methods. Finite difference methods are widely favored due to their simplicity, efficiency, and ease of implementation. For wave equations in anisotropic media, high-order finite difference schemes (such as regular grids, staggered grids, or rotated staggered grids) can be directly used to discretize the spatial partial derivatives, achieving high spatial simulation accuracy. For the discretization of time derivatives, second-order difference schemes are currently the main method used. These low-precision discretization schemes are prone to strong temporal numerical dispersion and even numerical instability. With the development of high-order finite difference methods, spatial high-order difference discretization has been widely used. However, existing time high-order difference schemes are primarily designed for wave equations in isotropic media. Due to their simpler form and fewer parameters involved, these numerical solutions cannot be directly extended to anisotropic media. The Lax-Wendroff time discretization scheme can be extended to anisotropic media, but its computational efficiency is relatively low. Other high-precision spectral methods, due to their involvement with Fourier transforms, struggle to maintain computational efficiency, especially for three-dimensional wave equations. By introducing paraxial mesh nodes into axial nodes, an improved difference template is constructed that achieves high-order accuracy in both the time and spatial domains for numerical simulations of isotropic wave equations.

[0004] Overall, 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 includes a time term high-order difference coefficient, and the space discrete function includes 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] Solving the simplified spatiotemporal dispersion relation based on a target optimization algorithm to obtain optimized high-order differential coefficients of the time term and optimized high-order differential coefficients of the space term;

[0012] extrapolating and solving the three-dimensional VTI medium acoustic wave equation based on the optimized high-order differential coefficient of the time term, the optimized high-order differential coefficient of the space term, the time discrete function, and the space discrete function;

[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), (-m,-n,-l), (-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), (-m,-n,-l), (-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, 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 above-described three-dimensional anisotropic medium time-space domain seismic forward modeling methods.

[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, comprising 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 difference template, which can effectively improve the time and space simulation accuracy; then, a simplified processing is performed and the solution is obtained through a target optimization algorithm to obtain 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 following briefly introduces the drawings required for use in the embodiments. 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 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 three-dimensional anisotropic media in time and space domains 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 1 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 N-Method (12,4) method 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 N-Method (12,4) method 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 difference method and the improved method in a three-dimensional uniform model at 0.25 s according to an embodiment of the present application; FIG5(a) is a wavefield snapshot calculated by the traditional difference method C-Method (12, 2) and Δt = 1.25 s for the three-dimensional uniform model at 0.25 s, FIG5(b) is a wavefield snapshot calculated by the traditional difference method C-Method (12, 2) and Δt = 0.625 s for the three-dimensional uniform model at 0.25 s, and FIG5(c) is a wavefield snapshot calculated by the method N-Method (12, 4) of the present application and Δt = 1.25 s for the three-dimensional uniform model at 0.25 s.

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

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

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

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

[0035] Anisotropy is widely present in the earth's medium. 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 technologies play an important role in seismic data processing. 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 wave equation of three-dimensional VTI (transverse isotropy with avertical axis of symmetry, i.e. VTI, vertical symmetry axis transverse isotropy) media, which can effectively improve the numerical accuracy in anisotropic wave field simulation.

[0036] In order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the present application is further described in detail below with reference to 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, integrated on the server 104, or 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, the server 104 performs an improved high-order discrete format derivation on the three-dimensional uniform VTI medium model to be processed, simplifies the derivation of the dispersion relationship and the calculation of the differential coefficients, and thus 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, specifically, it can be executed by a computer device such as a terminal or a server alone, or it can be executed 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, including 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 anisotropic earth medium model in the field of seismic exploration. In a three-dimensional VTI medium, the second-order wave equation characterized by the pressure component, i.e., 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 using the finite difference method to numerically solve the three-dimensional VTI medium acoustic wave equation, the traditional solution usually uses a high-order difference format to solve the spatial derivative. Solve the problem and apply the second-order difference scheme to the time derivative However, relevant numerical experiments have shown that when the time sampling interval is large, the traditional second-order time difference format is prone to cause strong time dispersion and even propagation instability. In addition, existing high-order time 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, this application provides a time-space high-order finite difference numerical solution in step 202 below 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 a time term high-order difference coefficient, and the space discrete function contains a space term high-order difference coefficient.

[0045] Specifically, relevant research results show that adding additional non-axial grid nodes to the traditional differential template 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), (-m,-n,-l), (-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), (-m,-n,-l), (-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, 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, based on the characteristics of plane waves, the present application expresses the wave field components in the three-dimensional VTI medium acoustic wave equation as follows:

[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 simplifies the time-discrete function and the space-discrete function based on the three-dimensional VTI medium acoustic wave equation to obtain a simplified time-space domain dispersion relation. Specifically, the equations corresponding to the preset improved time high-order difference operator and the preset improved space high-order difference operator are substituted into the three-dimensional VTI medium acoustic wave equation. After simplification, the original time-space domain dispersion relation for the coupling between the anisotropy parameter and the spatial wave number is obtained, 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] Note 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 we can get:

[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 with 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 combining G = 2cos(ωΔt)-2, we can get the simplified time-space domain dispersion relationship:

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

[0067] Where ω is the angular frequency.

[0068] Step 204 : Solve the simplified spatiotemporal dispersion relation based on the target optimization algorithm to obtain the optimized high-order differential coefficients of the time term and the optimized high-order differential coefficients of the space term.

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

[0070] (1) Taylor series is used to expand the trigonometric functions at both ends of the simplified time-space domain dispersion relation. In this step, the simplified time-space domain dispersion relation needs to be first 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 difference coefficients, g m is the arrangement combination of difference coefficients a 0,0,0 , a m,0,0 , a m,n,0 and a m,n,l ; χ m (β,θ,φ) has the following specific form in different parameter ranges:

[0080] When m≤M,

[0081] When M 2 / 4), m1=m2,

[0082]

[0083] When M 2 / 4), 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] where m1, m2, m3, m4, m5 are all mathematical operators, and all have no explicit physical meaning.

[0092] Considering that the Taylor expansion difference coefficient can only maintain high accuracy in the low frequency / low wave number area, in order to obtain high simulation accuracy in the entire wave number / frequency range, a target optimization algorithm is used to solve it here.

[0093] (2) Taking the minimum error at both ends of the preset function equation as the target, 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 based on 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 optimized high-order differential coefficients of the time term and the optimized high-order differential coefficients of the space term; specifically, using the least squares optimization method to solve the optimization objective function to obtain the optimized high-order differential coefficients of the time term and the optimized high-order differential coefficients of the 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 : extrapolating and solving the three-dimensional VTI medium acoustic wave equation 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.

[0101] In addition, 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] Here, κ 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. The other 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 suffers from significant 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 achieve higher numerical accuracy. Furthermore, in the optimized difference scheme, numerical accuracy can be further improved as the time operator N increases.

[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 adopt a smaller time step to ensure 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 adopt a larger time step to obtain high-precision simulation results, effectively balancing calculation accuracy and efficiency.

[0106] The application also uses a complex model to test the accuracy of different methods. FIG. 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 x 2500m x 1800m. The source is a Ricker wavelet, which is placed in the shallow layer of the model to generate vibration. The length of the spatial domain difference operator is 12. A mixed absorption boundary condition with a thickness of 10 layers is used to suppress artificial truncation boundary reflections. FIG. 8 shows the instantaneous wave field snapshots calculated by the traditional difference method and the method of the application using different time steps for the 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 difference method will exhibit unstable propagation. Therefore, in this example, the time step of the traditional difference method is 0.75ms, and the time step of the method of the application is 1.0ms. It can be observed Figure 8(a)-Figure 8(b) It can be seen that, in order to ensure the stability of wave field extrapolation, the traditional difference method needs to use a smaller time step, which inevitably increases the calculation time. In contrast, the method of the application considers the propagation characteristics of waves in the time domain and the spatial domain, and can effectively improve the numerical simulation accuracy while ensuring the calculation efficiency. Therefore, overall, the improved time high-order-space high-order finite difference numerical solution method in the application can provide an effective wave field continuation tool for subsequent anisotropic medium high-precision seismic processing.

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

[0108] In an example embodiment, a computer device is provided, which can be a server or a terminal. The computer device comprises a processor, a memory, an input / output interface (I / O) and a communication interface. The processor, the memory and the input / output interface are connected through a system bus, and the communication interface is connected to the system bus through the input / output interface. The processor of the computer device is configured to provide computing and control capabilities. The memory of the computer device comprises 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 running the operating system and the computer program in the non-volatile storage medium. The database of the computer device is configured to store a three-dimensional VTI medium acoustic wave equation. The input / output interface of the computer device is configured to exchange information between the processor and external devices. The communication interface of the computer device is configured to communicate with external terminals through a network connection. The computer program is executed by the processor to implement a three-dimensional anisotropic medium space-time domain seismic forward method.

[0109] In an example embodiment, a computer device is provided, which comprises a memory and a processor. The memory stores a computer program, and the processor executes the computer program to implement the steps in the above method embodiments.

[0110] In an example embodiment, a computer readable storage medium is provided, which stores a computer program. The computer program is executed by a processor to implement the steps in the above method embodiments.

[0111] In an example embodiment, a computer program product is provided, which comprises a computer program. The computer program is executed by a processor to implement the steps in the above method embodiments.

[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 for analysis, stored data, displayed data, etc.) involved in the present application are all information and data authorized by the user or authorized by all parties, and the collection, use and processing of related data need to comply with relevant regulations.

[0113] Those skilled in the art will understand that all or part of the processes in the above-mentioned embodiment methods can be implemented 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 memory, database or other media used in the embodiments provided in this application may include at least one of non-volatile and volatile memory. Non-volatile memory may 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 may 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 databases involved in the various embodiments provided herein may include at least one of a relational database and a non-relational database. Non-relational databases may include, but are not limited to, distributed databases based on blockchains. The processors involved in the various embodiments provided herein may include, but are not limited to, general-purpose processors, central processing units, graphics processing units, digital signal processors, programmable logic units, data processing logic units based on quantum computing, and the like.

[0115] The technical features of the above embodiments can be combined arbitrarily. In order 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 document uses specific examples to illustrate the principles and implementation methods of this application. The description of the above examples is only intended to help understand the method and core concept of this application. At the same time, for those skilled in the art, based on the concept of this application, there may 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 by: The three-dimensional anisotropic medium time-space domain seismic forward modeling method includes: 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 includes a time term high-order difference coefficient, and the space discrete function includes 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; Solving the simplified spatiotemporal dispersion relation based on a target optimization algorithm to obtain optimized high-order differential coefficients of the time term and optimized high-order differential coefficients of the space term; extrapolating and solving the three-dimensional VTI medium acoustic wave equation based on the optimized high-order differential coefficient of the time term, the optimized high-order differential coefficient of the space term, the time discrete function, and the space discrete function; 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), (-m,-n,-l), (-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), (-m,-n,-l), (-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, 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, 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 time-space domain 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, 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: Taylor series is used to expand the trigonometric functions at both ends of the simplified time-space domain dispersion relationship. 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 A linear combination 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 The permutation and combination of θ is the plane wave propagation angle, φ is the plane wave azimuth, where the plane wave is obtained by decomposing the pressure wave or 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, Constructing an optimization objective function with the goal of minimizing the error between both ends of the preset function equation; 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: Where E is the error value and b is the upper limit 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 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 according to any one of claims 1 to 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 according to any one of claims 1 to 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 according to any one of claims 1 to 6 is implemented.

Images