Seismic wave forward modeling numerical simulation method, device and system and storage medium

Abstract

The invention relates to the technical field of numerical simulation, and discloses a seismic wave forward numerical simulation method, device and system and a storage medium, and the method comprises the steps: determining a first objective function based on a traditional regular grid finite difference format three-dimensional scalar wave equation; based on the first objective function, establishing a second objective function of a three-dimensional effective vector wave number range; a regularization parameter is introduced into the second objective function, an optimal solution is determined by solving a regularization problem, and an optimized finite difference coefficient is obtained; and substituting the optimized finite difference coefficient into the second objective function to obtain the maximum value of the effective vector wave number range. According to the method, the obtained three-dimensional space is established on the basis of the minimum norm, the domain second-order partial derivative rule grid finite difference coefficient is obtained, and error accumulation in numerical simulation of the deep target can be better reduced.

Description

Technical Field

[0001] The present invention relates to the field of numerical simulation technology, and in particular to a seismic wave forward numerical simulation method, device, system and storage medium. Background Art

[0002] Seismic wave forward modeling is divided into two categories: physical simulation and numerical simulation. Physical simulation requires large initial investment and is difficult to make materials, but the results are very realistic and reliable; numerical simulation is simple, efficient and low-cost. At present, the numerical simulation methods of seismic waves can be roughly divided into: integral equation method, ray tracing method and wave equation method. The wave equation method gives full consideration to the geometric and dynamic characteristics of seismic waves, and its application range is wider and its prospects are broader. At present, 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 modeling and reverse time migration technology will become a very important and meaningful topic.

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

[0004] (1) The numerical dispersion caused by the finite difference method is unavoidable due to its gridding discretization of the medium. When using low-order difference approximation, coarse grids with large step sizes, or high-frequency sources, serious numerical dispersion will occur, thus affecting the accuracy of wave field extension. 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 differential 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 they cannot adapt to large time sampling intervals and numerical simulation of seismic waves in deep high-speed media. The computational 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 the migration noise, efficiently and accurately preserving the source wave field, and improving the calculation speed of wave field extension in 3D RTM are also the basis for achieving high-precision and high-efficiency 3D RTM.

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

[0008] The optimization algorithm is used to obtain new differential coefficients for forward numerical simulation of wave equations, that is, the differential coefficients are optimized before wave field extension. The differential coefficients and differential order remain unchanged during the entire process of wave field extension, which can effectively suppress spatial dispersion without increasing the amount of calculation. According to the type of objective function for optimizing differential coefficients, the optimization methods can be roughly divided into three categories: First, based on the maximum norm (L ∞ ) optimization method can obtain the widest effective wavenumber range, but it will bring a large accumulation error. This is because L ∞ The L2 norm only focuses on the maximum approximation error, but ignores the error distribution within the required wavenumber range. Second, the constant coefficient based on the L2 norm uses an optimization method, that is, reducing the sum of squares of the approximation error within the effective wavenumber range. Third, compared with the L2 and L ∞ Compared with the optimization method based on the L1 norm, the finite difference coefficients provided by the optimization algorithm based on the L1 norm have the lowest approximation error in the low and medium wavenumber range; therefore, the highest numerical simulation accuracy can be achieved even in long-duration simulations.

[0009] According to the number of spatial directions involved in the optimization, constant coefficient optimization methods can be 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 approximate errors of other propagation angles except the propagation angle considered. The second type is the multi-directional optimization method, which simultaneously optimizes the finite difference coefficients in all spatial directions, considers the differential approximation errors of all wave field propagation angles in space, and can effectively reduce the dispersion anisotropy while satisfying the differential 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] In order to solve the above problems, the present invention provides a seismic wave forward numerical simulation method, device, system and storage medium to meet the requirements of three-dimensional high-precision forward numerical simulation, establish a three-dimensional objective function for obtaining the finite difference coefficients of the second-order partial derivative regular grid in the spatial domain, and obtain the difference coefficients by solving the ADMM algorithm. Under the condition of an absolute error of 0.5E-5 allowable error, the three-dimensional multi-directional optimization only sacrifices a small effective wavenumber coverage range, and can obtain a higher difference approximation accuracy for medium and low wavenumbers, and can also weaken the anisotropy of spatial dispersion errors.

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

[0013] According to a first embodiment of the present invention, a method for numerical simulation of seismic wave forward modeling is provided, the method comprising:

[0014] Determine the first objective function based on the three-dimensional scalar wave equation in the traditional regular grid finite difference format;

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

[0016] Introducing a regularization parameter into the second objective function, determining an optimal solution by solving the regularization problem, and obtaining optimized finite difference coefficients;

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

[0018] Furthermore, the determining of the first objective function based on the three-dimensional scalar wave equation in the conventional regular grid finite difference format comprises:

[0019] The three-dimensional scalar wave equation in the traditional regular grid finite difference format is expressed as:

[0020]

[0021] 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;

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

[0023]

[0024] 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 mth 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 next m wave fields in the x direction;

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

[0026]

[0027] Where k represents the wave number;

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

[0029]

[0030] Where E(k) represents the dispersion error.

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

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

[0033]

[0034] 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; is the Laplace operator, expressed as:

[0035]

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

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

[0038]

[0039] Where P0 represents the plane wave amplitude, i represents the imaginary unit, and θ represents the wave number The angle with the z-axis, Indicates wave number The azimuth of ω represents the angular frequency;

[0040] Substituting formula (7) into formula (5) yields:

[0041]

[0042] 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

[0043]

[0044] Where E represents the sum of absolute errors;

[0045] In the selected amplitude spectrum, the dominant frequency fm The frequency corresponding to the peak value 1 / 64 is the cutoff frequency fc When , the corresponding cutoff wave number is expressed as:

[0046]

[0047] In the formula, k c represents the cutoff wave number, v min Indicates the minimum value of velocity in the model;

[0048] In formula (9), k, θ, After discretization, they are expressed as k i =i*k max / I、θ j =j*π / J, I , J, 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:

[0049]

[0050] Where E(i,j,g) represents the dispersion error in three-dimensional space;

[0051] Express formula (11) in matrix form:

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

[0053] In the formula, ||·||1 represents the L1 norm; A is an L×N 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:

[0054]

[0055] In the formula, n represents the nth order, x l,1 express L×1 The matrix of order, x l,2 express L×2 The matrix of order, x l,3 express L×3 Matrix of order.

[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 unit matrix;

[0060] The optimal solution is found by solving the regularized problem, namely:

[0061]

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

[0063]

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

[0065] Further, substituting the optimized finite difference coefficients into the second objective function to obtain the maximum value of the effective vector wave number range includes:

[0066] 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 <η, 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.

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

[0068] 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;

[0069] 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;

[0070] 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;

[0071] The range maximum value calculation unit is configured to substitute the optimized finite difference coefficient into the second objective function to obtain the maximum value of the effective vector wave number range.

[0072] According to a third technical solution of the present invention, a seismic wave forward numerical simulation system is provided, the system comprising: a memory for storing a computer program; and 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, a non-transitory computer-readable storage medium storing instructions is provided. When the instructions are executed by a processor, the method described above is executed.

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

[0075] The present invention establishes a three-dimensional space to obtain the domain second-order partial derivative regular grid finite difference coefficient based on the minimum norm, which is more conducive to reducing the error accumulation in the numerical simulation of deep targets. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

[0078] Figure 3 A flow chart of a seismic wave forward numerical simulation method according to an embodiment of the present invention is shown.

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

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

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

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

[0083] Figure 8 The three-dimensional uniform model Rb single point acceptance seismic record according to an embodiment of the present invention and the residual diagram with the analytical solution (a: t: 2.3s-2.5s; b: t: 3.8s-4.0s) are shown.

[0084] Fig. 9 A schematic diagram of the high and steep part of a three-dimensional overthrust model according to an embodiment of the present invention is shown.

[0085] Fig.10 The single-point acceptance seismic record of the three-dimensional uniform model Ra according to an embodiment of the present invention and the residual diagram with the analytical solution are shown (a: t: 2.3s-2.5s; b: t: 3.8s-4.0s).

[0086] Fig.11 The three-dimensional uniform model Rb single-point acceptance seismic record according to an embodiment of the present invention and the residual diagram with the analytical solution are shown (a: t: 2.3s-2.5s; b: t: 3.8s-4.0s).

[0087] Fig.12 A structural diagram of a seismic wave forward numerical simulation device according to an embodiment of the present invention is shown. DETAILED DESCRIPTION

[0088] The following describes the embodiments of the present invention by specific examples, and those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that the following embodiments and features in the embodiments can be combined with each other without conflict.

[0089] The specific implementation of the present invention is further described in detail below in conjunction with the drawings and examples.

[0090] In order to meet the requirements of three-dimensional high-precision forward numerical simulation, an embodiment of the present invention provides a method for forward numerical simulation of seismic waves, that is, to establish a three-dimensional objective function for obtaining the finite difference coefficients of the second-order partial derivative regular grid in the spatial domain and to obtain the differential coefficients by solving the ADMM algorithm. Under the condition of an absolute error of 0.5E-5 allowable error, the three-dimensional multi-directional optimization only sacrifices a small effective wave number coverage range, and can obtain a higher differential approximation accuracy for medium and low wave numbers, while also reducing the anisotropy of the spatial dispersion error. It is known from the numerical simulation of three-dimensional uniform medium models and three-dimensional complex models that as the extension time of the wave field increases, the error accumulation phenomenon will become more serious, and the control of error accumulation is particularly important at this time. In order to be superior to unidirectional optimization in suppressing spatial numerical dispersion and reducing numerical dispersion anisotropy, a method is proposed, that is, to establish a three-dimensional space to obtain the second-order partial derivative regular grid finite difference coefficients of the domain 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 They are the dispersion error curves after single-direction optimization (ADMM-1D) and three-dimensional multi-directional optimization (ADMM-3D) under the maximum allowable error threshold η=0.5E-5, Δx=Δy=Δz=20m. From the dispersion curves of different differential orders, it can be seen that a wider effective wavenumber coverage range can be obtained by increasing the differential order.

[0091] Specifically, see Figure 3 , is a flow chart of a seismic wave forward numerical simulation method, the seismic wave forward numerical simulation method includes the following steps S100 to S400, which are described in detail as follows:

[0092] S100, determining a first objective function based on a three-dimensional scalar wave equation in a traditional regular grid finite difference format.

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

[0094]

[0095] Where P(x,y,z,t) is the pressure wave field, s(x,y,z,t) is the source function, v is the wave velocity, 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-order derivative is defined as:

[0097]

[0098] 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 mth 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.

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

[0100]

[0101] Where k represents the wave number;

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

[0103]

[0104] Where E(k) represents the dispersion error.

[0105] S200. Based on the first objective function, a second objective function of a three-dimensional effective vector wave number range is established.

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

[0107]

[0108] 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, P0,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] is the Laplace operator, expressed as:

[0110]

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

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

[0113]

[0114] Where P0 represents the plane wave amplitude, i represents the imaginary unit, and θ represents the wave number The angle with the z-axis, Indicates wave number The azimuth of ω represents the angular frequency;

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

[0116]

[0117] 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

[0118]

[0119] Where E represents the sum of absolute errors.

[0120] In the selected amplitude spectrum, the dominant frequency fm The frequency corresponding to the peak value 1 / 64 is the cutoff frequency fc When , the corresponding cutoff wave number is expressed as:

[0121]

[0122] In the formula, k c represents the cutoff wave number, v min Represents the minimum value of velocity in the model.

[0123] In formula (9), k, θ, After discretization, they are expressed as k i =i*k max / I、θ j =j*π / J, I , J, 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:

[0124]

[0125] Where 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, 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:

[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, introducing a regularization parameter into the second objective function, determining an optimal solution by solving the regularization problem, and obtaining optimized finite difference coefficients.

[0134] In order 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] Where α represents the regularization parameter and D represents the unit matrix.

[0137] The optimal solution is found by solving the regularized problem, namely:

[0138] The optimal solution is found by solving the regularized problem, namely:

[0139]

[0140] Equation (15) is a nonlinear problem. Currently, many methods can be used to solve it. However, in order to obtain constant coefficients efficiently and accurately, the alternating direction multiplier method is used to convert the original problem into a constraint problem:

[0141]

[0142] d=Ac+b, 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:

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

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

[0148] Step 3,

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

[0150] S400, substituting the optimized finite difference coefficients into the second objective function to obtain the maximum value of the effective vector wave number range.

[0151] When the differential order in the three directions and the maximum allowable error threshold are determined, in order to find the optimal solution, a smaller k can be set first. max , will produce an ε max ; When ε max <η, slightly increase k max Open and repeat the above process until the cycle ends. Table 1 shows the differential coefficients obtained by solving equation (11) using the ADMM algorithm with different differential orders.

[0152] when When and θ are not zero, it can be seen that the dispersion error anisotropy of single-direction optimization in the low and medium wavenumber domain is more obvious, and the overall error is larger; Figure 4As shown in the figure, the anisotropy of dispersion error is also quite obvious when the step sizes in the three directions are different (Δx = Δy = 20m, Δz = 10m, and the differential coefficients are shown in Table 2); however, the three-dimensional multi-directional optimization only needs to sacrifice a small effective wavenumber coverage range in exchange for a higher differential approximation accuracy in the medium and low wavenumber domains, which can effectively reduce the anisotropy of dispersion error and reduce the accumulation of errors, thereby improving the accuracy of wavefield 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] In order to verify the reliability of the method proposed in the present invention, forward numerical simulations were carried out on a three-dimensional uniform model and a three-dimensional overthrust model, and accuracy analysis and comparison were performed.

[0159] (1) Homogeneous medium model

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

[0161] from Figure 5 , Figure 6The error comparison of the mid-wave field snapshot shows that in the numerical simulation results of the two different optimization methods, the error control of the differential pair of three-dimensional multi-directional optimization is better than that of the single-directional optimization, which is consistent with the results of the dispersion relationship analysis; the wave field error value at the wave field snapshot t = 3.2s is larger than that at t = 0.4s, indicating that with the increase of wave field propagation time, the error accumulation is more serious, and the three-dimensional multi-directional optimization differential coefficient has a better control on the error accumulation. From the record analysis of the Ra point and the Rb point ( Figure 7 , Figure 8 ) can also draw this conclusion.

[0162] (2) Complex Model

[0163] Complex model This embodiment selects the high and steep part of the three-dimensional overthrust model such as Fig. 9 As shown in Figure 1, the model size is 4km×4km×4km, its P-wave velocity varies from 2200m / s to 6000m / s, the grid spacing is Δx=20m, Δy=20m, Δy=10m, the time step is 1ms, the Ricker wavelet with a main frequency of 25Hz is used, and the source is placed at the center of the model (2000m, 2000m, 2000m). As with the uniform model, without adding boundary conditions, the total sampling time is set to 4s. Fig.10 , Fig.11 The differences between the partial receiving records and analytical solutions of the single point Ra (3800m, 2000m, 2000m) and Rb (3000m, 3000m, 3200m) of the numerical simulation of the finite difference acoustic wave equation were obtained for different optimization methods. It can also be seen that the error control of the differential pair of three-dimensional multi-directional optimization is better than that of one-dimensional optimization, and the suppression effect of propagation error accumulation is better, which can further improve the accuracy of numerical simulation when simulating deep and complex geological models.

[0164] The embodiment of the present invention also provides a seismic wave forward numerical simulation device, such as Fig.12 As shown, the device comprises:

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

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

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

[0168] The range maximum value calculation unit 124 is configured to substitute the optimized finite difference coefficient into the second objective function to obtain the maximum value of the effective vector wave number range.

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

[0170] An embodiment of the present invention further provides a seismic wave forward numerical simulation system, the system comprising: a memory for storing a computer program; and 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, and when the instructions are executed by a processor, the method described in any of the above embodiments is executed.

[0172] The above implementation modes are only used to illustrate the present invention, but not to limit the present invention. Ordinary technicians in the relevant technical field can 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 should 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; The optimized finite difference coefficients are substituted into the second objective function to obtain the maximum value of the effective vector wave number range.

2. The seismic wave forward numerical simulation method according to claim 1, characterized in that: 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 order of 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 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.

3. The seismic wave forward numerical simulation method according to claim 2, 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, and z directions with different discrete step lengths, and 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; 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: Where P0 represents the plane wave amplitude, i represents the imaginary unit, and θ represents the wave number The angle with the z-axis, Indicates wave number The azimuth of ω 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 fm The frequency corresponding to the peak value 1 / 64 is the cutoff frequency fc 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 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 express L×1 The matrix of order, x l,2 express L×2 The matrix of order, x l,3 express L×3 Matrix of order.

4. The seismic wave forward numerical simulation method according to claim 3, 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.

5. The seismic wave forward numerical simulation method according to claim 4, 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 <η, 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.

6. 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; The range maximum value calculation unit is configured to substitute the optimized finite difference coefficient into the second objective function to obtain the maximum value of the effective vector wave number range.

7. 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 5. 8 . A non-transitory computer-readable storage medium storing instructions, which, when executed by a processor, executes the method according to claim 1 .

Images