A numerical forward simulation method, device, system and storage medium for seismic waves

Abstract

The present invention relates to the technical field of numerical simulation, and discloses a seismic wave forward numerical simulation method, device, system and storage medium. The method includes: determining a first objective function based on the three-dimensional scalar wave equation of the traditional regular grid finite difference format; establishing a second objective function within the three-dimensional effective vector wave number range based on the first objective function; introducing a regularization parameter into the second objective function, and determining an optimal solution by solving the regularization problem to obtain optimized finite difference coefficients; substituting the optimized finite difference coefficients into the second objective function to obtain the maximum value within the effective vector wave number range. The present invention establishes the second-order partial derivative regular grid finite difference coefficients in the three-dimensional space solution domain obtained based on the minimum norm, which is more conducive to reducing the error accumulation in the numerical simulation of deep targets.

Description

Technical Field

[0001] The present invention relates to the technical field of numerical simulation, and particularly relates to a forward seismic wave numerical simulation method, device, system and storage medium. Background Technique

[0002] The forward simulation of seismic waves is divided into two categories: physical simulation and numerical simulation. Physical simulation has large upfront investment and difficult material production, but the result authenticity and reliability are very high; numerical simulation is simple, efficient and low-cost. Currently, seismic wave numerical simulation methods can be roughly divided into: integral equation method, ray tracing method and wave equation method. The wave equation method fully considers the geometric and dynamic characteristics of seismic waves, and its application range is larger and the prospect is broader. Currently, most scholars use two-dimensional media for testing when studying the numerical simulation and reverse time migration of seismic waves. However, the limitations of two-dimensional simulation make it difficult to simulate the seismic response characteristics of actual underground three-dimensional geological bodies. Therefore, with the development of computer software and hardware, and the increase in the complexity of seismic exploration, the research on three-dimensional high-precision wave equation forward and reverse time migration technologies will become a very important and meaningful topic.

[0003] The finite difference numerical simulation technology is widely used in the forward numerical simulation of wave equations due to its high calculation efficiency and simple operation, and has become the basis for reverse time migration and even full waveform inversion. However, when performing forward numerical simulation and reverse time migration of three-dimensional finite difference (FD) method wave equations, the following aspects are worthy of attention:

[0004] (1) Due to the grid discretization of the medium by the finite difference method, the numerical dispersion brought about is inevitable. Serious numerical dispersion will occur when using low-order difference approximation, coarse grids with large step sizes or high-frequency seismic sources, which will affect the accuracy of wave field continuation. Therefore, for the numerical simulation of three-dimensional seismic waves, it is very important to control the error accumulation of wave field propagation.

[0005] (2) The difference coefficients solved by conventional methods cannot obtain a large maximum Courant-Friedrichs-Lewy (CFL) number that satisfies the stability condition, that is, the stability factors are not high, which means that it cannot adapt to the seismic wave numerical simulation of large time sampling intervals and deep high-velocity media. The calculation efficiency of the forward wave field and the memory overhead required for wave field storage in reverse time migration will be directly affected.

[0006] (3) Effectively suppressing migration noise, efficiently and accurately storing the source wave field and improving the calculation speed of wave field continuation in three-dimensional reverse time migration are also the basis for achieving high-precision and high-efficiency three-dimensional reverse time migration.

[0007] The finite-difference (FD) method is popular in the numerical simulation of wave propagation due to its simple implementation and low computational cost. However, due to the approximation of the derivative operator by the finite-difference operator, the FD method has inherent numerical errors. Therefore, the numerical phase velocity of wave propagation varies with frequency and propagation angle. The time-dispersion transformation proposed by Koene can completely remove the time dispersion generated under the time-step condition far exceeding the stability limit. Therefore, this paper mainly focuses on suppressing spatial numerical dispersion and weakening the numerical dispersion error anisotropy caused by the spatial finite-difference operator.

[0008] Use an optimization algorithm to obtain new difference coefficients for the forward numerical simulation of the wave equation, that is, optimize the difference coefficients before wave-field extrapolation, and the difference coefficients and difference orders remain unchanged throughout the wave-field extrapolation process. In this way, both spatial dispersion can be effectively suppressed and the computational amount will not increase. According to the type of objective function for optimizing the difference coefficients, the optimization methods can be roughly divided into three categories: First, the optimization method based on the maximum norm (L ∞ ) can obtain the widest effective wave-number range, but it will bring large cumulative errors because the L ∞ norm only focuses on the maximum approximation error and ignores the error distribution within the required wave-number range. Second, the constant-coefficient optimization method based on the L2 norm, that is, reducing the sum of the squares of the approximation errors within the effective wave-number range. Third, compared with the optimization methods based on the L2 and L ∞ norms, the optimization algorithm based on the L1 norm provides the finite-difference coefficients with the lowest approximation errors in the medium and low wave-number ranges; therefore, even in long-duration simulations, the highest numerical simulation accuracy can be obtained.

[0009] According to the number of spatial directions involved in the optimization, the constant-coefficient optimization method can be further divided into two categories. The most common type is the single-direction optimization method, which only considers the finite-difference coefficients in a single spatial direction and cannot take into account the dispersion anisotropy because it ignores the approximation errors at other propagation angles except the considered propagation angle. The second type is the multi-direction optimization method, which simultaneously optimizes the finite-difference coefficients in all spatial directions, considers the difference approximation errors at all wave-field propagation angles in space, and can effectively weaken the dispersion anisotropy under the condition of meeting the difference approximation accuracy.

[0010] However, the existing forward numerical simulation methods cannot meet the requirements of three-dimensional high-precision forward numerical simulation. Summary of the Invention

[0011] To solve the above problems, the present invention provides a seismic wave forward numerical simulation method, apparatus, system, and storage medium, which can meet the requirements of three-dimensional high-precision forward numerical simulation, establish an objective function for obtaining the second-order partial derivative regular grid finite difference coefficients in the three-dimensional spatial domain, and solve the obtained difference coefficients through the ADMM algorithm. Under the allowable error condition of an absolute error of 0.5E-5, the three-dimensional multi-directional optimization can sacrifice only a tiny effective wave number coverage range to obtain a higher difference approximation accuracy for medium and low wave numbers, and at the same time can also weaken the anisotropic phenomenon of spatial dispersion error.

[0012] To achieve the above object, the technical solution adopted by the present invention is as follows:

[0013] According to the first aspect of the present invention, a seismic wave forward numerical simulation method is provided, and the method includes:

[0014] Based on the three-dimensional scalar wave equation of the traditional regular grid finite difference format, determine the first objective function;

[0015] Based on the first objective function, establish a second objective function for the three-dimensional effective vector wave number range;

[0016] Introduce a regularization parameter into the second objective function, and determine the optimal solution by solving the regularization problem to obtain the optimized finite difference coefficients;

[0017] Substitute the optimized finite difference coefficients into the second objective function to obtain the maximum value of the effective vector wave number range.

[0018] Further, the determining the first objective function based on the three-dimensional scalar wave equation of the traditional regular grid finite difference format includes:

[0019] Express the three-dimensional scalar wave equation of the traditional regular grid finite difference format as:

[0020]

[0021] In the formula, v represents the wave speed, P represents the pressure wave field, t represents time, x represents the x direction, y represents the y direction, z represents the z direction, and s(t) represents the source function;

[0022] In the x direction, the 2M-order finite difference staggered grid of the spatial second derivative is defined as:

[0023]

[0024] In the formula, Δx represents the discrete step length in the x direction, c m (m = 1, 2,..., M) represents the finite difference coefficient, m represents the m-th difference coefficient, M represents the last difference coefficient, 2M represents the difference approximation order, c0 represents the initial finite difference coefficient, P0,0,0 represents the initial wave field, P -m,0,0 represents the first m wave fields in the x direction, P m,0,0 represents the last m wave fields in the x direction;

[0025] Based on the one-dimensional plane wave theory, the wave number is calculated by the following formula:

[0026]

[0027] In the formula, k represents the wave number;

[0028] With the goal of obtaining the finite difference coefficients, a first objective function is established, and the first objective function is expressed as:

[0029]

[0030] In the formula, E(k) represents the dispersion error.

[0031] Furthermore, based on the first objective function, a second objective function for the three-dimensional effective vector wave number range is established, including:

[0032] Using different discrete step sizes in the x, y, and z directions, a new spatial partial derivative is established, expressed as:

[0033]

[0034] In the formula, respectively represent the finite difference coefficients of the second-order partial derivative difference approximation in the x, y, and z directions, M x 、M y 、M z represent the lengths of the finite difference operators in the x, y, and z directions; △y and △z respectively represent the discrete step sizes in the y and z directions; P 0,m,0 represents the last m wave fields in the y direction, P 0,-m,0 represents the first m wave fields in the y direction, P 0,0,m represents the last m wave fields in the z direction, P 0,0,-m represents the first m wave fields in the z direction; is the Laplace operator, expressed as:

[0035]

[0036] In the formula, P(x, y, z, t) represents the pressure wave field;

[0037] The expression of the plane wave theory in three-dimensional space is determined as follows:

[0038]

[0039] where \(P_0\) represents the plane wave amplitude, \(i\) represents the imaginary unit, \(\theta\) represents the angle between the wave number \(k\) and the \(z\)-axis, \(\varphi\) represents the azimuth angle of the wave number \(k\), and \(\omega\) represents the angular frequency;

[0040] Substituting Equation (7) into Equation (5) gives:

[0041]

[0042] According to Equation (8), the differential coefficients are optimized by minimizing the sum of the absolute errors within the given effective wave number \(k\) max range, and the sum of the absolute errors is expressed as

[0043]

[0044] where \(E\) represents the sum of the absolute errors;

[0045] When selecting the main frequency \(f\) in the amplitude spectrum m and the corresponding frequency of the peak 1 / 64 is the cut-off frequency \(f_c\) c , the corresponding cut-off wave number is expressed as:

[0046]

[0047] where \(k_c\) c represents the cut-off wave number, and \(v_{min}\) min represents the minimum value of the velocity in the model;

[0048] After discretizing \(k\), \(\theta\), and in Equation (9), they are respectively expressed as \(k_i\) i = \(i\cdot k / I\), \(\theta_j\) max = \(j\cdot\pi / J\), and j \(\varphi_g\) = \(g\cdot2\pi / G\), where \(I\), \(J\), and \(G\) respectively represent the discretization degrees of \(k\), \(\theta\), and \(\varphi\), and the first objective function is rewritten as a second objective function in the three-dimensional effective vector wave number range, expressed as: where \(E(i,j,g)\) represents the dispersion error in the three-dimensional space;

[0049]

[0050]

[0051] Equation (11) is expressed in matrix form as:

[0052] \(E(c)=\|Ac + b\|_1\ (12)\)

[0053] where \(\|\cdot\|_1\) represents the \(L_1\) norm; \(A\) is an \(L\times N\) matrix, \(L = I\cdot J\cdot G\); \(b\) is an \(L\)-dimensional vector; and \(c\) is the optimized \(N\) finite difference coefficients, with the specific form as follows:

[0054] ​

[0055] In the formula, n represents the nth order, and x l,1 represents an L×1 order matrix, and x l,2 represents an L×2 order matrix, and x l,3 represents an L×3 order matrix.

[0056] Furthermore, a regularization parameter is introduced into the second objective function, and the optimal solution is determined by solving the regularization problem to obtain the optimized finite difference coefficients, including:

[0057] Introducing the regularization parameter into the second objective function, we get:

[0058]

[0059] In the formula, ψ(c) represents, F(c) represents, α represents the regularization parameter, and D represents the identity matrix;

[0060] Finding the optimal solution by solving the regularization problem, that is:

[0061]

[0062] d = Ac + b, and the corresponding augmented Lagrangian function is expressed as:

[0063]

[0064] In the formula, l(c, d) represents the augmented Lagrangian function, η represents the maximum allowable error, and u represents the soft threshold.

[0065] Furthermore, substituting the optimized finite difference coefficients into the second objective function to obtain the maximum value of the effective vector wavenumber range, including:

[0066] When determining the difference orders in three directions and the maximum allowable error threshold conditions, set an initial k max , based on the generated ε max , when ε max < η, increase k max And repeat the above process until the loop ends, completing the solution of the second objective function to obtain the maximum value of the effective vector wavenumber range.

[0067] According to the second technical solution of the present invention, a seismic wave forward numerical simulation device is provided, and the device includes:

[0068] A first objective function determination unit configured to determine a first objective function based on the three-dimensional scalar wave equation in the traditional regular grid finite difference format;

[0069] A second objective function determination unit, configured to establish a second objective function for a three-dimensional effective vector wavenumber range based on the first objective function;

[0070] A finite difference coefficient optimization unit, configured to introduce a regularization parameter into the second objective function, and determine an optimal solution by solving a regularization problem to obtain optimized finite difference coefficients;

[0071] A range maximum value calculation unit, configured to substitute the optimized finite difference coefficients into the second objective function to obtain the maximum value of the effective vector wavenumber range.

[0072] According to a third technical solution of the present invention, there is provided a seismic wave forward numerical simulation system, the system including: a memory for storing a computer program; a processor for executing the computer program to implement the method as described above.

[0073] According to a fourth technical solution of the present invention, there is provided a non-transitory computer-readable storage medium storing instructions, which when executed by a processor, execute the method as described above.

[0074] The present invention has at least the following beneficial effects:

[0075] The present invention establishes a second-order partial derivative regular grid finite difference coefficient for a three-dimensional space acquisition domain obtained based on the minimum norm, which is more conducive to reducing error accumulation in numerical simulation of deep targets. Description of the Drawings

[0076] Figure 1 Shows a one-dimensional optimized dispersion curve diagram according to an embodiment of the present invention (η = 0.5E-5, Δx = Δy = Δz = 20m).

[0077] Figure 2 Shows a three-dimensional optimized dispersion curve diagram according to an embodiment of the present invention (η = 0.5E-5, Δx = Δy = Δz = 20m).

[0078] Figure 3 Shows a flowchart of a seismic wave forward numerical simulation method according to an embodiment of the present invention.

[0079] Figure 4 Shows a comparison diagram of dispersion curves of different optimization methods according to an embodiment of the present invention (Mx = My = 8, Mz = 4η = 0.5E-5, Δx = Δy = 20, Δz = 10m).

[0080] Figure 5 Shows a wave field snapshot diagram at t = 0.4s according to an embodiment of the present invention (a: analytical solution; b: residual between ADMM-1D and analytical solution; c: residual between ADMM-3D and analytical solution).

[0081] Figure 6 Shows the t = 3.2 s wave field snapshot according to an embodiment of the invention (a: analytical solution; b: residual between ADMM-1D and analytical solution; c: residual between ADMM-3D and analytical solution).

[0082] Figure 7 Shows the single-point received seismic record of the three-dimensional homogeneous model Ra according to an embodiment of the present invention and the residual map with the analytical solution (a: t: 2.3 s to 2.5 s; b: t: 3.8 s to 4.0 s).

[0083] Figure 8 Shows the single-point received seismic record of the three-dimensional homogeneous model Rb according to an embodiment of the present invention and the residual map with the analytical solution (a: t: 2.3 s to 2.5 s; b: t: 3.8 s to 4.0 s)

[0084] Figure 9 Shows the schematic diagram of the high-steep part of the three-dimensional reverse thrust nappe model according to an embodiment of the present invention.

[0085] Figure 10 Shows the single-point received seismic record of the three-dimensional homogeneous model Ra according to an embodiment of the present invention and the residual map with the analytical solution (a: t: 2.3 s to 2.5 s; b: t: 3.8 s to 4.0 s).

[0086] Figure 11 Shows the single-point received seismic record of the three-dimensional homogeneous model Rb according to an embodiment of the present invention and the residual map with the analytical solution (a: t: 2.3 s to 2.5 s; b: t: 3.8 s to 4.0 s).

[0087] Figure 12 Shows the structural diagram of a seismic wave forward numerical simulation device according to an embodiment of the present invention. Detailed implementation manners

[0088] The following uses specific specific examples to illustrate the implementation manners of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific implementation manners. Various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, without conflict, the following embodiments and the features in the embodiments can be combined with each other.

[0089] The following further describes in detail the specific implementation manners of the present invention in conjunction with the drawings and embodiments.

[0090] To meet the requirements of three-dimensional high-precision forward numerical simulation, an embodiment of the present invention provides a seismic wave forward numerical simulation method, that is, to establish an objective function for obtaining the second-order partial derivative regular grid finite difference coefficients in the three-dimensional spatial domain and solve the obtained difference coefficients through the ADMM algorithm. Under the allowable error condition of an absolute error of 0.5E-5, the three-dimensional multi-directional optimization can obtain higher difference approximation accuracy for medium and low wave numbers while sacrificing only a small effective wave number coverage range, and can also weaken the anisotropy phenomenon of spatial dispersion error. From the numerical simulations of the three-dimensional homogeneous medium model and the three-dimensional complex model, it is known that as the wave field continuation time increases, the error accumulation phenomenon will become more serious, and at this time, the control of error accumulation becomes particularly important. In order to be superior to the one-way optimization in suppressing the spatial numerical dispersion and weakening the numerical dispersion anisotropy, a method is proposed, that is, to establish the second-order partial derivative regular grid finite difference coefficients in the three-dimensional spatial domain obtained based on the minimum norm, which is more conducive to reducing the error accumulation in the numerical simulation of deep targets. Figure 1 and Figure 2 are the dispersion error curves after the single-direction optimization (ADMM-1D) and the three-dimensional multi-directional optimization (ADMM-3D) respectively under the maximum allowable error threshold η = 0.5E-5, Δx = Δy = Δz = 20m. It can be seen from the dispersion curves of different difference orders that a wider effective wave number coverage range can be obtained by increasing the difference order.

[0091] Specifically, please refer to Figure 3 , which is a flowchart of a seismic wave forward numerical simulation method. The seismic wave forward numerical simulation method includes the following steps S100 to step S400, which are introduced in detail as follows:

[0092] S100. Determine the first objective function based on the three-dimensional scalar wave equation of the traditional regular grid finite difference format.

[0093] In this embodiment, according to the three-dimensional scalar wave equation of the traditional regular grid finite difference format, it can be expressed as:

[0094]

[0095] In the formula, P(x, y, z, t) is the pressure wave field, s(x, y, z, t) is the source function, v is the wave speed, P represents the pressure wave field, t represents time, x represents the x direction, y represents the y direction, z represents the z direction, and s(t) represents the source function.

[0096] Taking x as an example, the 2M-order finite difference staggered grid of the spatial second derivative is defined as:

[0097]

[0098] In the formula, Δx represents the grid step size in the x direction, c m(m = 1, 2, …, M) represents the finite difference coefficient, m represents the m-th difference coefficient, M represents the last difference coefficient, 2M represents the order of the difference approximation, c0 represents the initial finite difference coefficient, P 0,0,0 represents the initial wave field, P -m,0,0 represents the first m wave fields in the x direction, P m,0,0 represents the last m wave fields in the x direction.

[0099] Since P(x, y, z, t) is a scalar, according to the one-dimensional plane wave theory P(x) = P0e ikx It can be obtained that:

[0100]

[0101] In the formula, k represents the wave number;

[0102] In order to obtain the difference coefficient, the first objective function that can be established is as follows:

[0103]

[0104] In the formula, E(k) represents the dispersion error.

[0105] S200. Based on the first objective function, establish the second objective function for the three-dimensional effective vector wave number range.

[0106] In this embodiment, since the wave number k (k x , k y , k z ) is a vector, |k| = k, therefore, the difference coefficient obtained by using formula (4) will make the spatial dispersion error anisotropy in the non-axis direction more serious and cannot adapt to the discrete mode under non-cubic grid conditions. In order to adapt to the undulating ground surface and the high-steep structure model, different discretization step lengths in the x, y, and z directions can be adopted, and the new spatial partial derivative can be expressed as:

[0107]

[0108] In the formula, respectively represent the finite difference coefficients of the second-order partial derivative difference approximation in the x, y, and z directions, M x , M y , M z represent the lengths of the finite difference operators in the x, y, and z directions; △y and △z respectively represent the discretization step lengths in the y and z directions; P 0,m,0 represents the last m wave fields in the y direction, P 0,-m,0 represents the first m wave fields in the y direction, P 0,0,m represents the last m wave fields in the z direction, P 0,0,-m represents the first m wave fields in the z direction.

[0109] The Laplace operator is denoted as:

[0110]

[0111] In the formula, P(x, y, z, t) represents the pressure wave field.

[0112] The expression for determining the plane wave theory in three-dimensional space is as follows:

[0113]

[0114] In the formula, P0 represents the plane wave amplitude, i represents the imaginary unit, θ represents the angle between the wave number k and the z-axis, represents the azimuth angle of the wave number k, and ω represents the angular frequency;

[0115] Substituting formula (7) into formula (5) gives:

[0116]

[0117] According to formula (8), within the given effective wave number k max range, the difference coefficients are optimized by minimizing the total absolute error, and the total absolute error is denoted as

[0118]

[0119] In the formula, E represents the total absolute error.

[0120] When selecting the main frequency f m in the amplitude spectrum and the peak 1 / 64 corresponds to the cut-off frequency f c , the corresponding cut-off wave number is denoted as:

[0121]

[0122] In the formula, k c represents the cut-off wave number, and v min represents the minimum value of the velocity in the model.

[0123] After discretizing k, θ, in formula (9), they are respectively denoted as k i = i * k max / I, θ j = j * π / J, I , J, and G respectively represent the discretization degrees of k, θ, The first objective function is rewritten as a second objective function in the three-dimensional effective vector wave number range, denoted as:

[0124]

[0125] In the formula, E(i, j, g) represents the dispersion error in three-dimensional space.

[0126] For the convenience of calculation, formula (11) is expressed in the form of a matrix:

[0127] E(c) = ||Αc + b||1 (12)

[0128] In the formula, ||·||1 represents the L1 norm; A is an L×N matrix, where L = I*J*G; b is an L-dimensional vector; c is the optimized N finite difference coefficients, and the specific form is as follows:

[0129]

[0130] In the formula, n represents the nth order, x l,1 represents an L×1 matrix, x l,2 represents an L×2 matrix, x l,3 represents an L×3 matrix.

[0131] x L×3 is an L×3 matrix, which can be obtained by the following calculation:

[0132]

[0133] S300. Introduce the regularization parameter into the second objective function, and determine the optimal solution by solving the regularization problem to obtain the optimized finite difference coefficients.

[0134] To overcome the numerical oscillation that may be caused by directly optimizing the constant coefficient c (Wang, 2016), the regularization technique (Tikhonov, 2013) can be introduced into the objective function to obtain:

[0135]

[0136] In the formula, α represents the regularization parameter, and D represents the identity matrix.

[0137] Find the optimal solution by solving the regularization problem, that is:

[0138] Find the optimal solution by solving the regularization problem, that is:

[0139]

[0140] Equation (15) is a non-linear problem. Currently, many methods can be used to solve it. However, in order to obtain the constant coefficient efficiently and accurately, the alternating direction multiplier method is used to transform the original problem into a constrained problem:

[0141]

[0142] d = Ac + b, and the corresponding augmented Lagrangian function can be written as:

[0143]

[0144] where l(c, d) represents the augmented Lagrangian function, η represents the maximum allowable error, and u represents the soft threshold.

[0145] Exemplarily, the update steps for regularizing and solving the objective function based on the above method are as follows:

[0146] Step 1, Input: A, b, the maximum number of iterations K, the maximum allowable error η;

[0147] Step 2, Set the optimization parameters: c 0 = 0, u 0 = 0, ρ, k = 0;

[0148] Step 3,

[0149] Step 4, Output: the optimized finite difference coefficients.

[0150] S400. Substitute the optimized finite difference coefficients into the second objective function to obtain the maximum value of the effective vector wavenumber range.

[0151] When the difference orders in three directions and the maximum allowable error threshold conditions are determined, in order to find the optimal solution, a relatively small k can be set first max , which will generate an ε max ; when ε max < η, slightly increase k max and repeat the above process until the loop ends. Table 1 shows the difference coefficients obtained by solving Equation (11) using the ADMM algorithm for different difference orders.

[0152] When and θ are both non-zero, it can be seen that the dispersion error anisotropy in the medium and low wavenumber domain is relatively obvious for single-direction optimization, and the overall error is relatively large; as Figure 4 shown, when the step sizes in three directions are different, the dispersion error anisotropy is also relatively obvious (Δx = Δy = 20m, Δz = 10m, and the difference coefficients are shown in Table 2); however, three-dimensional multi-directional optimization only needs to sacrifice a tiny effective wavenumber coverage range to obtain a higher difference approximation accuracy in the medium and low wavenumber domain, can effectively reduce the dispersion error anisotropy, reduce the accumulation of errors, and thus improve the accuracy of wave field simulation.

[0153] Table 1 Finite Difference Coefficients Solved by ADMM-3D

[0154]

[0155]

[0156] Table 2 Differential coefficients of different optimization methods (Mx = 8, My = 8, Mz = 4)

[0157]

[0158] To verify the reliability of the method proposed in the present invention, forward numerical simulations are respectively carried out on a three-dimensional homogeneous model and a three-dimensional inverse thrust nappe model below, and accuracy analysis and comparison are carried out.

[0159] (1) Homogeneous medium model

[0160] In the homogeneous medium model, the longitudinal wave velocity is 2000 m / s, the model size is 4 km × 4 km × 4 km, the grid spacing is Δx = 20 m, Δy = 20 m, Δy = 10 m, the time step is 1 ms, a Ricker wavelet with a dominant frequency of 25 Hz is used, and the source is placed at the center position (2000 m, 2000 m, 2000 m) of the model. To obtain the propagation record of the acoustic wavelength duration, the total sampling time is set to 4 s without applying boundary conditions. The differences between the wavefield snapshots of the numerical simulation of the finite-difference acoustic wave equation obtained by different optimization methods and the analytical solution (Green's function solution) are compared as Figure 5 、 Figure 6 shown Figure 7 、 Figure 8 are respectively the differences between the single-point partial reception records at Ra (3760 m, 2000 m, 2000 m) and Rb (2960 m, 2960 m, 3120 m) and the analytical solution.

[0161] From Figure 5 、 Figure 6 the error comparison of the wavefield snapshots, it can be seen that among the numerical simulation results of the two different optimization methods, the differential of the three-dimensional multi-directional optimization is better than that of the single-directional optimization in terms of error control, which is consistent with the results of the dispersion relation analysis; the wavefield error value at the moment of t = 3.2 s in the wavefield snapshot is larger than that at the moment of t = 0.4 s, indicating that as the wavefield propagation time increases, the error accumulation is more serious, and the three-dimensional multi-directional optimization differential coefficient has better control over the error accumulation. This conclusion can also be drawn from the analysis of the reception records at points Ra and Rb( Figure 7 、 Figure 8 ).

[0162] (2) Complex model

[0163] In this complex model, the high-steep part of the three-dimensional inverse thrust nappe model is selected in this embodiment as Figure 9As shown, the model size is 4 km × 4 km × 4 km, the P-wave velocity varies from 2200 m / s to 6000 m / s, the grid spacing is Δx = 20 m, Δy = 20 m, Δy = 10 m, the time step is 1 ms, a Ricker wavelet with a dominant frequency of 25 Hz is used, and the source is placed at the center of the model (2000 m, 2000 m, 2000 m). Similar to the homogeneous model, without applying boundary conditions, the total sampling time is set to 4 s. Figure 10 , Figure 11 They are the differences between the partial receiving records of single points Ra(3800 m, 2000 m, 2000 m) and Rb(3000 m, 3000 m, 3200 m) obtained by numerical simulation of the finite-difference acoustic wave equation using different optimization methods and the analytical solution, respectively. It can also be seen that the three-dimensional multi-directional optimized difference is superior to the one-dimensional optimization in error control and has a better suppression effect on the accumulation of propagation errors. When simulating a deep and complex geological model, it can further improve the numerical simulation accuracy.

[0164] The embodiment of the present invention also provides a seismic wave forward numerical simulation device, as Figure 12 shown, the device includes:

[0165] A first objective function determination unit 121, configured to determine a first objective function based on the three-dimensional scalar wave equation in the traditional regular grid finite-difference format;

[0166] A second objective function determination unit 122, configured to establish a second objective function within the three-dimensional effective vector wave number range based on the first objective function;

[0167] A finite-difference coefficient optimization unit 123, configured to introduce a regularization parameter into the second objective function, solve the regularization problem to determine the optimal solution, and obtain the optimized finite-difference coefficient;

[0168] A range maximum value calculation unit 124, configured to substitute the optimized finite-difference coefficient into the second objective function to obtain the maximum value within the effective vector wave number range.

[0169] It should be noted that the device described in this embodiment and the method described above belong to the same technical concept, have the same technical principle, and can achieve the same beneficial effects, so it will not be elaborated here.

[0170] The embodiment of the present invention also provides a seismic wave forward numerical simulation system, the system includes: a memory for storing a computer program; a processor for executing the computer program to implement the method described in any of the above embodiments.

[0171] An embodiment of the present invention further provides a non-transitory computer-readable storage medium storing instructions, which, when executed by a processor, perform the method described in any one of the above embodiments.

[0172] The above embodiments are only used to illustrate the present invention and are not intended to limit the present invention. Those of ordinary skill in the relevant technical field can also make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, all equivalent technical solutions also belong to the scope of the present invention. The patent protection scope of the present invention shall be defined by the claims.

Claims

1. A seismic wave forward numerical simulation method, characterized in that: The method comprises: Determine the first objective function based on the three-dimensional scalar wave equation in the traditional regular grid finite difference format; Based on the first objective function, establishing a second objective function of a three-dimensional effective vector wave number range; Introducing a regularization parameter into the second objective function, determining an optimal solution by solving the regularization problem, and obtaining optimized finite difference coefficients; Substituting the optimized finite difference coefficients into the second objective function to obtain the maximum value of the effective vector wave number range; The method of determining the first objective function based on the three-dimensional scalar wave equation in the conventional regular grid finite difference format comprises: The three-dimensional scalar wave equation in the traditional regular grid finite difference format is expressed as: In the formula, v represents wave velocity, P represents pressure wave field, t represents time, x represents x direction, y represents y direction, z represents z direction, and s(t) represents source function; In the x-direction, the 2M-order finite difference staggered grid of the spatial second-order derivative is defined as: Where Δx represents the discrete step length in the x direction, c m ,m=1,2,…,M represents the finite difference coefficient, m represents the first difference coefficient, M represents the last difference coefficient, 2M represents the difference approximation order, c0 represents the initial finite difference coefficient, P 0,0,0 represents the initial wave field, P -m,0,0 represents the first m wave fields in the x direction, P m,0,0 represents the next m wave fields in the x direction; Based on one-dimensional plane wave theory, the wave number is calculated by the following formula: Where k represents the wave number; With the goal of obtaining the finite difference coefficient, a first objective function is established, and the first objective function E(k) is expressed as: Where E(k) represents the dispersion error.

2. The seismic wave forward numerical simulation method according to claim 1, characterized in that: The step of establishing a second objective function of a three-dimensional effective vector wave number range based on the first objective function includes: use x、y、z Different discrete step lengths in three directions establish new spatial partial derivatives, expressed as: In the formula, They represent the finite difference coefficients of the second-order partial derivative approximation in the x, y, and z directions, respectively. x 、M y 、M z represents the length of the finite difference operator in the x, y, and z directions; △y and △z represent the discrete step lengths in the y and z directions respectively; P 0,m,0 represents the last m wave fields in the y direction, P 0,-m,0 represents the first m wave fields in the y direction, P 0,0,m represents the last m wave fields in the z direction, P 0,0,-m Represents the first m wave fields in the z direction; ▽ 2 P is the Laplace operator, expressed as: Where P(x,y,z,t) represents the pressure wave field; The expression that determines the plane wave theory in three-dimensional space is as follows: In the formula, P0 represents the plane wave amplitude, i represents the imaginary unit, θ represents the angle between the wave number k and the z-axis, represents the azimuth of wave number k, ω represents the angular frequency; Substituting formula (7) into formula (5) yields: According to formula (8), at a given effective wave number k max Minimize the sum of absolute errors within the range to optimize the differential coefficients. The sum of absolute errors is expressed as Where E represents the sum of absolute errors; In the selected amplitude spectrum, the dominant frequency f m The frequency corresponding to the peak value 1 / 64 is the cutoff frequency f c When , the corresponding cutoff wave number is expressed as: In the formula, k c represents the cutoff wave number, v min Indicates the minimum value of velocity in the model; In formula (9), k, θ, After discretization, they are expressed as k i =i*k max / I、θ j =j*π / J, I, J, and G represent k, θ, The discreteness of the first objective function is rewritten into the second objective function of the three-dimensional effective vector wave number range, which is expressed as: Where E(i,j,g) represents the dispersion error in three-dimensional space; Express formula (11) in matrix form: E(c)=||Αc+b||1 (12) In the formula, ||·||1 represents the L1 norm; A is an l×n-order matrix, l=I*J*G; b is an l-order vector; c is the optimized n finite difference coefficients, and the specific form is as follows: In the formula, n represents the nth order, x l,1 represents an l×1 matrix, x l,2 represents an l×2-order matrix, x l,3 Represents an l×3 matrix.

3. The seismic wave forward numerical simulation method according to claim 2, characterized in that: The regularization parameter is introduced into the second objective function, and the optimal solution is determined by solving the regularization problem to obtain the optimized finite difference coefficients, including: Introducing the regularization parameter into the second objective function, we get: In the formula, α represents the regularization parameter, and D represents the unit matrix; The optimal solution is found by solving the regularized problem, namely: d=Ac+b, the corresponding augmented Lagrangian function is expressed as: Where l(c, d) represents the augmented Lagrangian function, η represents the maximum allowable error, and u represents the soft threshold.

4. The seismic wave forward numerical simulation method according to claim 3, characterized in that: Substituting the optimized finite difference coefficients into the second objective function to obtain the maximum value of the effective vector wave number range includes: When the differential order and the maximum allowable error threshold in the three directions are determined, an initial k is set. max , based on the generated ε max , when ε max When <η, increase k max The above process is repeated until the loop ends, the second objective function is solved, and the maximum value of the effective vector wave number range is obtained.

5. A seismic wave forward numerical simulation device, characterized in that: The device comprises: A first objective function determination unit is configured to determine a first objective function based on a three-dimensional scalar wave equation in a conventional regular grid finite difference format; A second objective function determining unit is configured to establish a second objective function of a three-dimensional effective vector wave number range based on the first objective function; A finite difference coefficient optimization unit is configured to introduce a regularization parameter into the second objective function, determine an optimal solution by solving a regularization problem, and obtain optimized finite difference coefficients; a range maximum value calculation unit, configured to substitute the optimized finite difference coefficient into the second objective function to obtain a maximum value of the effective vector wave number range; The method of determining the first objective function based on the three-dimensional scalar wave equation in the conventional regular grid finite difference format comprises: The three-dimensional scalar wave equation in the traditional regular grid finite difference format is expressed as: In the formula, v represents wave velocity, P represents pressure wave field, t represents time, x represents x direction, y represents y direction, z represents z direction, and s(t) represents source function; In the x-direction, the 2M-order finite difference staggered grid of the spatial second-order derivative is defined as: Where Δx represents the discrete step length in the x direction, c m ,m=1,2,…,M represents the finite difference coefficient, m represents the first difference coefficient, M represents the last difference coefficient, 2M represents the difference approximation order, c0 represents the initial finite difference coefficient, P 0,0,0 represents the initial wave field, P -m,0,0 represents the first m wave fields in the x direction, P m,0,0 represents the next m wave fields in the x direction; Based on one-dimensional plane wave theory, the wave number is calculated by the following formula: Where k represents the wave number; With the goal of obtaining the finite difference coefficient, a first objective function is established, and the first objective function E(k) is expressed as: Where E(k) represents the dispersion error.

6. A seismic wave forward numerical simulation system, characterized in that: The system comprises: Memory for storing computer programs; A processor, configured to execute the computer program to implement the method according to any one of claims 1 to 4. 7 . A non-transitory computer-readable storage medium storing instructions, which, when executed by a processor, executes the method according to claim 1 .

Images