Reducing the low-wavenumber dispersion error by building the Lagrange dual problem with a powerful local restriction

Abstract

1. Introduction

The finite-difference (FD) method (Alterman & Karal 1968; Kelly et al. 1976) is an efficient approach to simulate seismic wave propagation in a discrete medium because of its low computational memory consumption and simple numerical implementation. Nevertheless, the FD scheme causes a severe numerical error when using a coarse spatial grid. Additionally, some seismic exploration techniques, such as least-squares migration imaging (LeBras & Clayton 1988; Nemeth et al.1999) and full waveform inversion (Tarantola 1984, 1987; Yong et al.2019), require a significant amount of computational cost. Therefore, geophysicists desire an efficient and accurate numerical simulation method.

The algorithms to compute the FD coefficients can be categorized as the Taylor-series expansion (TE) method and optimization algorithms. The TE method easily obtains the FD coefficients by using Taylor-series expansion of dispersion relation in particular domain, but not able to handle the numerical dispersion well. Optimizing FD coefficients is a more efficient approach to improve computational efficiency and suppress numerical dispersion, compared with FD coefficients of the TE method.

The TE methods can be briefly divided into two types: TE in the space domain (Dablain 1986; Fornberg 1998) and TE in the time-space domain (Liu & Sen 2009). Under the same length of FD coefficients as other methods, although the TE methods present a smaller numerical dispersion in the low-wavenumber region, they possess larger dispersion in the high-wavenumber region (Holberg 1987; Etgen and O’Brien 2007).

For reducing this significant error in the high-wavenumber region, Holberg (1987) first proposed reducing the maximal error of the group velocity within a limited band realized by optimizing FD coefficients. Liu (2013, 2014) then minimized the error of the dispersion relations in the space domain using the least-square (LS) method, which can obtain the global solution without iteration. The cost function is the crux to the optimal FD coefficients for the optimization algorithm, directly affecting in the ability to suppress the dispersion error. Liu (2013) proved that the cost function based on the relative error of dispersion relation results in a lower dispersion error within the effective bandwidth than that based on the absolute error. Therefore, we set the relative error as the cost function in this paper. The dispersion error contains the spatial dispersion error and the temporal dispersion errors; this paper’s dispersion error is the spatial dispersion error.

For the optimization problem, the solution of the cost function based on the |${L}_\infty $| norm has a wider effective bandwidth than that of the |${L}_2$| norm cost function (Zhang & Yao 2013a,b). The algorithms involved in the solving process of the |${L}_\infty $| norm cost function include the simulated annealing method (SA, Zhang & Yao 2013a,b), Remez exchange algorithm (REA, Kosloff et al.2010; An 2015; Yang et al.2016, 2017; He et al. 2019) and LS (Holberg 1987). Using the optimal FD coefficients to high accurately simulate the broad seismic wave propagation accurately was proposed by Zhang & Yao 2013a,b). Due to the improper initial parameters and restrictions, the SA algorithm still has a narrower effective bandwidth, compared with REA. Appling the REA to obtain the optimal FD coefficients was an efficient way for solving the |${L}_\infty $| norm cost function in the past decade. Kosloff et al. (2010) used the REA (McClellan & parks 1972; Soubaras 1996) to solve an |${L}_\infty $| norm cost function for obtaining the implicit filter coefficients. Yang et al. (2016) proposed a new implicit FD scheme that applies the sampling approximation method to approach the dispersion relation for obtaining FD coefficients with broad effective bandwidth. Also, the sampling approximation method by Yang et al. (2016) is similar to the first calculation in the Remez exchange method. An (2015) introduced a Remez-like exchange method to obtain the uniform approximation scheme of the dispersion relation for all wavenumbers in the frequency domain. Regarding the staggered grid FD scheme, and given the production error by the narrow effective bandwidth of the conventional FD coefficients numerically, applying REA for computing optimal FD coefficients can improve this phenomenon (Yang et al.2017). He et al. (2018) modified the near-zero wavenumber condition by regarding the |${c}_1$| (2M order FD scheme contain M + 1 FD coefficients, |${c}_0,{c}_1,\ \ldots \ $|⁠, |${c}_1$| is the second FD coefficient) as an unknown, which leads to the equal-ripple amplitude error property; thus, the new method has the widest effective bandwidth in the current optimization strategies.

The FD coefficients should obey two criteria, the first one is that the optimal FD coefficients should apply to the largest possible bandwidth within the acceptable error, the second one is that it should possess the smallest possible dispersion error over a wide bandwidth. For the first criterion, the |${L}_\infty $| norm optimization problem solved by REA has the widest bandwidth within the given error tolerance. Regarding the second one, Miao & Zhang (2020) verify that the |${L}_1$| norm cost function has a more robust capability to suppress dispersion error than employing |${L}_2$| and |${L}_\infty $| norms. However, the ADMM solved |${L}_1$| norm cost function only considers the overall wavenumbers in the effective bandwidth (Miao & Zhang 2020), which lacks local constraints on the low wavenumbers. We contrapose this issue and propose our optimal method with a powerful restriction.

This paper analyses the dispersion curves of the three optimization algorithms and summarizes the relationship between the wavenumbers of zero points and dispersion error in the low-wavenumber region. Based on this relationship, the sampling approximation method reduces the dispersion error, however, this method cannot adaptively adjust the effective bandwidth within a given error tolerance. To overcome this disadvantage, the establishment of the Lagrange dual problem is an efficient approach, which contains the |${L}_1$| norm cost function and the modified wavenumber restricted condition. The |${L}_1$| norm cost function is the same as that of Miao & Zhang (2020). The modified wavenumbers restricted condition aims to intensify the suppression of error in appointed wavenumber region, built by the modified wavenumbers of the REA zero points. The physical meaning of this restricted condition is that the dispersion errors of the new method at these modified wavenumbers are equal to 0. We introduce the dual ascent method to solve the Lagrange dual problem. Through analysis and simulations, we further demonstrate the superiority of the proposed optimization strategy, compared with methods based on the |${L}_1$| (Miao & Zhang 2020), |${L}_2$| (Liu 2013) and |${L}_\infty $| (He et al. 2018) norm-based cost functions.

2. Theory

2.1. The dispersion relation of the 2D acoustic wave equation

The (2M)th-order FD scheme of the spatial second derivatives can be written as (Kelly et al.1976; Dablain 1986)

where M is half of the FD order; |$u_m^n$| is a scalar wave field, for the second-order temporal and 2M order spatial FD scheme; the m and n are the loop counters |${\rm{m}} \in [ {0,{\rm{\ M}}} ]$|⁠, |${\rm{n}} \in ( { - 1,0,1} )$|⁠, |${\rm{m}},{\rm{\ n}} \in {\rm{Z}}$|⁠; h is the spatial step and |${c}_m$| are the FD coefficients.

According to the plane wave theory (Liu & Sen 2009), we have

where x, and t are the independent variables in the spatial domain and temporal domains, respectively, k denotes the wavenumber, |${\rm{\omega }}$| denotes the angular frequency, |$i{\rm{\ }}$|represents an imaginary number and |${\rm{\tau }}$| is the temporal step. Substituting equation (2) into equation (1) provides the dispersion relation in the space domain

where |$\beta = kh$| and maximal value is π, according to the Nyquist’s theorem. When |$\beta = 0$|⁠, the relation between the |${c}_0$| and |${c}_m$| can be expressed as

Substituting equation (4) into (3), the dispersion relation can be rewritten as

Equation (5) is known as the dispersion relation.

2.2. A series of special wavenumbers

The cost function to optimize the FD coefficients is always constructed by the dispersion relation. The cost function is divided into two types. The first type is the relative error of dispersion relation, expressed as

and the second type is the absolute error of dispersion relation, expressed as

For an optimal FD coefficient obtained by ADMM, LS, SA or REA, we found there are M wavenumbers that hold equation (6) to equal zero and M + 1 wavenumbers hold equation (7) to equal zero, as shown in figure 1a and b. The β = 0 is a wavenumber perpetually met in equation (7) to equal 0, thus, to get rid of β = 0, there are M wavenumbers still making equation (7) equal 0. Actually, the two sets of M wavenumbers are identical. For example, there is a wavenumber making the equation (6) equal to 0, which can be given as |$0 = \frac{{2\sum\nolimits_{m = 1}^M {{c}_m[ {\cos ( {m{\beta }_1} ) - 1} ]} }}{{\beta _{^1}^2}} + 1$|⁠, which can be easily changed into |$0 = 2\sum\nolimits_{m = 1}^M {{c}_m[ {\cos ( {m{\beta }_1} ) - 1} ]} + \beta _{^1}^2$|⁠. So, we refer to these wavenumbers as special M wavenumbers defined as |$\beta \_zer{o}_i( {i = 1,2,3...,M} )$|⁠.

Figure 1.

Comparison of the relative and absolute error for the REA, LS, SA and ADMM, as provided by He et al. (2019), Liu (2013, 2014), Miao & Zhang (2020) and Zhang & Yao 2013a,b).

Open in new tabDownload slide

As presented in figure 2 and Table 1, when the special M wavenumbers are concentrated in a low-wavenumber region, the FD coefficients possess a small dispersion error. Based on this, we compute the modified special M wavenumbers using the following formula

where p is a variable that alters the special M wavenumbers to focus on the low-wavenumber region. Applying the modified special M wavenumbers obtained from equation (8) to the sampling approximation method can reduce the low-wavenumber dispersion error, and the linear system of the sampling approximation method can be expressed as

Figure 2.

Relative error curve and the position of special M wavenumber for three methods, according to He et al. (2019), Liu (2013, 2014) and Miao & Zhang (2020).

Open in new tabDownload slide

The entries of equation (9) are as follows

We assume the row rank of A is smaller than M, and two linearly dependent row vectors must exist. The two row vectors can be represented as

where the |${A}_P$| is the Pth row vector of A, |${A}_K$| is the Kth row vector of A, |$K\& P = 1,2...M$|⁠, |$T \in R$|⁠.

Expanding equation (11), we can obtain

where |${\beta }_i \in (0,\pi ]\& {\beta }_i \in Q$|⁠. To make equation (12) hold, all |$\cos ( {P \times {\beta }_i} )$| and |$\cos ( {K \times {\beta }_i} )$| should be equal to 1. Further, all |$P \times {\beta }_i$| and |$K \times {\beta }_i$| should be equal to |$2n\pi $|⁠. Then, there must exist M|${\beta }_i$| to make the formulation hold, and the formulation is represented as

We using the maximal j\! (M) to compute the |${\beta }_i$|⁠, however, we found only M/2 |${\beta }_i$| exist, and thus the equation (13) does not hold. Therefore, the assumption that the row rank of A is smaller than M does not hold. So, the row rank of A is equal to M, and since the A is a square matrix, the rank of A is M.

We have proved the rank of A is M, so the row vector |${A}_j( {j = 1,2,...M} )$| of A is linearly independent, and the |${A}_j$| is |$\{ {[ {\cos ( {1{\beta }_j} ) - 1} ],[ {\cos ( {2{\beta }_j} ) - 1} ],...,[ {\cos ( {M{\beta }_j} ) - 1} ]} \}.$|

The augmented matrix (A d) is combined by the A and d, where A is an M*M matrix and d is a column vector with M elements. So, the augmented matrix (A d) is an M*M + 1 matrix. The rank of the augmented matrix (A d) is equal to the number of non-zero rows of the augmented matrix.

If the row vector of the augmented matrix (A d) is linearly independent, the row vector of the augmented matrix (A d) cannot be reduced to the zero vector by the row transformation, the number of non-zero rows of the augmented matrix (A d) will be equal to M. The row vector of the augmented matrix (A d) can be expressed as |${B}_j$|⁠, the |${B}_j$| is |$\{ {[ {\cos ( {1{\beta }_j} ) - 1} ],[ {\cos ( {2{\beta }_j} ) - 1} ],...,[ {\cos ( {M{\beta }_j} ) - 1} ], - \frac{{\beta _j^2}}{2}} \}$|⁠. However, according to the theorem (He et al.2014), when vector |${B}_j$| is obtained by adding one element into linearly independent vector |${A}_j$|⁠, the vector |${B}_j$| also is linearly independent. Thus, the row vector of the augmented matrix (A d) is linearly independent; furthermore, the rank of the augmented matrix (A d) is equal to M.

The rank of A and the rank of the augmented matrix (A d) are equal to M. Hence, equation (10) has only one solution according to the theory of linear equations. In other words, a set of special M wavenumbers can only derive one set of optimal FD coefficients. As presented in figure 3 and Table 2, the dispersion error of the new optimal FD coefficients is less than the ADMM under the same error limitation. Nevertheless, the new approach to optimize FD coefficients cannot adaptively adjust the error limitation and effective bandwidth. This opposes LS (Liu 2013), ADMM (Miao & Zhang 2020) and REA (He et al. 2019), which adjust the effective bandwidth according to the given error limitation. The adaptive adjustment of the error limitation and effective bandwidth means that a method can be easily used for numerical modeling and contrast research. Introducing an optimization problem with the L₁ norm cost function and the restricted condition is a good way to realize the adaptive adjustment of the effective bandwidth and dispersion error. Thus, the dual Lagrange problem is the best choice.

Figure 3.

Comparison of the relative error curve of new FD coefficients and ADMM.

Open in new tabDownload slide

2.3. The dual Lagrange problem

The L₁ norm cost function in the dual Lagrange problem is based on the relative error over the given range |$[ {0,{\beta }_{max}} ]$|

The numerical calculation of equation (14) can be considered as a sum of integrating a series of small intervals, the length of small intervals is |${{\rm{\beta }}}_{max}/N$|⁠, where N = 1000. The cost function |${\rm{E}}\_{\rm{in}}{{\rm{t}}}_{relative}$| is then the same as the cost function in Miao & Zhang (2020) by ignoring the constant item |${{\rm{\beta }}}_{max}/N$|⁠, which can be expressed as

Rearranging equation (15) yields the linear system of equations:

where the |$\| \cdot {\|}_1$| denotes the |${L}_1$| model, c is the optimal FD coefficients, A is the coefficient matrix and b is the coefficient vector, whose entries are as follows

|${\rm{w}}$|here |${{\rm{A}}}_{i,j}( {i = 1,2,3,\ \ldots \ N,j = 1,2,3,\ \ldots \ ,M} )$| represents the element in the coefficient matrix A, whose location is the ith row and the jth column and |${{\rm{b}}}_i$| denotes to the ith element in the vector b. Equation (16) is the L₁ norm cost function.

The Lagrange dual problem can be written as

The |${R}_{relative}( {{\bf c}} ) = 0$| is constructed by a modified special M wavenumber. The process of modifying the special M wavenumbers can be expressed as

where the |${\rm{\beta }}\_{\rm{zer}}{{\rm{o}}}_i$| is the special M wavenumber of optimal FD coefficients obtained by REA, the |$p = \beta \_zer{o}_1$|⁠. How to obtain these |${\rm{\beta }}\_{\rm{zer}}{{\rm{o}}}_i$| is shown in detail in the next section, and why this value is chosen is discussed in Appendix A. The modified wavenumbers restricted condition in the Lagrange dual problem can be written as

The physical meaning of equation (21) is that the dispersion error at the modified special M wavenumbers equals 0. Equation (21) can also be expressed in the form of the L₁ norm:

where |${{{\bf A}}}_1,{{\bf c}},{{{\bf b}}}_1$| are

|${\rm{w}}$|here |${{{\bf A}}}_{1_{i,j}}\!( {i = 1,2,3,\ \ldots \ N,j = 1,2,3,\ \ldots \ ,M} )$| represents the element in the coefficient matrix |${{{\bf A}}}_1$|⁠, whose location is the ith row and the jth column, and |${{\rm{b}}}_i$| denotes the ith element in the vector |${{{\bf b}}}_1$|⁠.

The cost function and the modified wavenumbers restricted condition combine as

where the |$\lambda $| is the Lagrange multiplier.

When the cost function is a convex function, the equality constraints are an affine function, the inequality constraints are a convex function and the optimization problem is a convex optimization problem. The norm function is a convex function according to Boyd’s theorem (Boyd & Vandenberghe 2004), and the affine function is the function constructed by the polynomial of the first order. Thus, the Lagrange dual problem (equations (18) and (19)) thus is a convex optimization problem. For a convex optimization problem without inequality constraint, this problem naturally holds a strong duality.

These features ensure that equation (24) can be solved by a dual ascent method (Boyd et al.2011) and further guarantee the uniqueness of the solution to the optimization problem. We apply the dual ascent to solve this optimization problem and obtain the optimal FD coefficients efficiently and accurately.

After several iterations of the two subproblems, the dual ascent method can solve the optimal problem. These subproblems contain updating c and updating the Lagrange multiplier |$\lambda $|⁠.

The first subproblem is updating c, where the new c can be expressed as

and |${A}_2 = [ {A:{\lambda }^k{A}_1} ]$| and |${b}_2 = {[ {{b}^T:{\lambda }^kb_1^T} ]}^T$|⁠. This is a standard LS problem, the minimum is obtained at the derivative vanish, so, the solution of this problem can be solved by the |$\frac{{\partial f}}{{\partial c}} = 0$|⁠. In addition, the optimal |${{{\bf c}}}^{k + 1}$| satisfies

The |${{{\bf c}}}^{k + 1}$| is obtained by the inverse of |${{\bf A}}_2^{{\bf T}}{{{\bf A}}}_2$|⁠. However, the solution of the inverse matrix requires a large computing resource. Therefore, the solution of equation (27) can be determined by introducing the conjugate-gradient method or Newton method.

Based on the Boyd et al. (2011) proposed dual ascent method principle, the second subproblem involves updating the Lagrange variable |$\lambda $|⁠, the new |$\lambda $| can be expressed as

where|$\ \lambda $| is the Lagrange variable and the |${\alpha }^k$| is the iteration direction. Here, we use the initial |$\lambda = 40$| and |${\alpha }^k$| = 10. The two steps of dual ascent and the closed solution of the new |$\lambda $| are mentioned in section 2.2 in Boyd et al.(2011). The initial value of c can be the optimal FD coefficients obtained by another method, such as the LS, REA, ADMM or TE. After a few iterations of equation step 1 and step 2, the error of the cost function converges to a minimum. As mentioned before, we have shown that the optimization problem has the only solution, and the solution will converge to the global solution regardless of the initial |$\lambda $| and the |${\alpha }^k$| values. However, the Lagrange multiplier λ and |${\alpha }^k$| should not be too small to avoid needing multiple iterations. According to our testing and development experience, the initial λ should exceed 40 and |${\alpha }^k$| should be >10 to ensure fewer than 100 iterations.

For the ADMM (Miao & Zhang 2020), three highly complex subproblems should be solved. However, the dual ascent method only computes two subproblems, which improve the method’s repeatability. Here we present this new and efficient method as an alternative to the geophysics community. Note that this method also can solve the L₁ cost function based on the absolute error of dispersion relation (Appendix A).

2.3. The Remez exchange algorithm

The REA is an approach that quickly obtains the best-approximation polynomial |$f( {\beta \_reme{z}_i} )$| of a given continuous function |$p( {\beta \_reme{z}_i} )$|⁠. The relationship between the |$f( {\beta \_reme{z}_i} )$| and |$p( {\beta \_reme{z}_i} )$| can be expressed as

Consider that the wavenumber domain dispersion relation, |$p( {\beta \_reme{z}_i} )$| is 1, and |$f( {\beta \_reme{z}_i} ) = \frac{{2\mathop \sum \nolimits_{m = 1}^M {c}_m[ {1 - cos( {m\beta \_reme{z}_i} )} ]}}{{\beta \_remez_i^2}}$|⁠. The optimal FD coefficients obtained by REA satisfy the equal-ripple property of the dispersion error magnitude, which possesses the widest bandwidth among the existing optimal FD coefficients (He et al. 2019) according to the Chebyshev criteria.

We then define the following error function by equation (3):

Combining equations (29) and (30), the equal-ripple error is represented as

where F is a parameter that changes as the iteration proceeds, |$\beta \_reme{z}_i$| is the set of wavenumbers, |$\beta \_reme{z}_i \in [ {0,k{h}_{max}} ]$| and |$k{h}_{max}$| is the maximum value of optimized wavenumber region. By rearranging equation (31), the matrix-vector equation computed the optimal FD coefficients can be expressed as

To obtain |$\beta \_zer{o}_i$| in the Lagrange dual problem, we derive the optimal FD coefficients of the REA. The existing research (He et al. 2019) determined the iterating process of REA. Once the REA is completed, the |$\beta \_zer{o}_i\ ( {i = 1,2,3,\ \ldots \ ,M} )$| in equation (17) can be obtained from the dispersion curve of REA FD coefficients.

3. Theoretical error analysis

The relative error of the dispersion relation in the wavenumber domain is written as

In this experiment analysis, the theoretical error is calculated by equation (33), which aims to demonstrate that the dispersion error of the proposed method is smaller than ADMM (L₁ norm cost function, Miao & Zhang 2020) within the effective bandwidth. The optimal FD coefficients of each approach are reported in Table 3. For the 16th-order FD coefficients, we calculate the relative error of three methods, e.g. the proposed method, the ADMM with the same effective bandwidth as the proposed method and the ADMM with the same error limitation as the proposed method. The corresponding results are illustrated in figure 4, which, under the same effective bandwidth or error limitation, demonstrates that the proposed method has a lower relative error in the black box than the ADMM. Hence, the solution of the Lagrange dual problem has a smaller numerical simulation error than ADMM. However, under the same given limitation, the effective bandwidth of the proposed method is slightly narrower than the ADMM, which is because of the exiting tradeoff between the effective bandwidth improvement and the dispersion error reduction in the field of optimal FD coefficients (Liu 2013; He et al.2018; Miao & Zhang 2020).

Figure 4.

Relative error of the 16th-order FD operators among the different optimization methods. The given error limitation is 0.0001.

Open in new tabDownload slide

4. Numerical simulations

We will solve the acoustic equation using a homogeneous example and the Marmousi model. The FD scheme with the 16th-order optimal FD coefficients discretizes the continuous acoustic equation, which is computed by the proposed method, the |${L}_\infty $| norm cost function (He et al. 2019), the |${L}_2$| norm cost function (Liu 2013) and the |${L}_1$| norm cost function (Miao & Zhang 2020) under the same given error tolerance (0.0001). Besides, the Ricker wavelet is a source-time function with a dominant frequency of 20 Hz. The velocity of the homogeneous square model is 1500 m s^-1, the grid system is 2000 × 2000, the spatial step is 5 m, and the temporal step is 0.0001 s and to avoid the accumulation of the temporal dispersion, the propagation time is 1, 2 and 3 s. We chose the snapshots and the residual snapshots as the result of this experiment, and to conveniently observe the residual error of each method, we extract the single trace at x = 5000 m for a demonstration. The size sequence of residual error in snapshots is the REA, least-square methods (LSM), ADMM and proposed method, as illustrated in figures 5 and 6, which indicate that the key factor to determining the residual error is the dispersion error of optimal FD coefficients. Furthermore, figures 5 and 6 illustrate that the enhancive dispersion error within the bandwidth will generate the numerical dispersion. Obviously, figures 5 and 6 highlight that the proposed method has a smaller dispersion error than the competitor methods regarding the homogeneous model simulation test. With the increased propagation time and distance, the residual error of each method is accumulated, conforming to the error rule in the FD scheme. Figure 6 infers that the proposed method has the lowest velocity of accumulated error, indicating the feasibility of the proposed method for the long-time seismic wave simulation.

Figure 5.

Combined snapshots and wavefield errors for four methods in the homogeneous model. The first column refers to snapshots at 1, 2 and 3 s, respectively. The second column is the residual error between the reference solution and snapshots calculated by four types of the optimal FD operator. The reference solution is computed by the 40th-order FD coefficients from the Taylor expansion method in the space domain.

Open in new tabDownload slide

Figure 6.

Traces from figure 5b, d and f, in which the traces with the x coordinate of 5000 m from are from (a) figure 5b; (b) figure 5d and (c) figure 5f. We highlight the specific peak error on the y-axis for all optimized methods.

Open in new tabDownload slide

To verify the proposed method’s applicability to seismic techniques involved various velocity models, e.g. the LS reverse time migration or full waveform inversion production, we design an experiment that uses the different velocity model simulations to check the advantage of the proposed method. Hence, we compare the ADMM and the proposed method, as the ADMM is superior to LSM and REA (figure 6). The velocities of the homogeneous models are 3000 and 5000 m s⁻¹, respectively. The model size is 2000 × 2000 m, the spatial step is 10 m, the temporal step is 0.0002 and the propagation distance of the seismic wave is 600 m. The dominant 20 Hz frequency Ricker wavelet is located at (1000 × 1000 m). We thus use the trace of the simulative result at 1000 m to reveal the precision of the ADMM and the proposed method, as depicted in figure 7. Figure 7 illustrates that the proposed method has a higher accurate modeling result than ADMM either in the 3000 or 5000 m s⁻¹ model simulation. This suggests that the proposed method is good at maintaining high precision in the various velocity model, and with the seismic technique with the iterative velocity model the proposed method does not require recomputing the FD coefficients.

Figure 7.

Trace comparison of the simulation result in the homogenous models with the different velocities for the 16th-order FD coefficients of the ADMM and proposed approach. The reference solution is computed by the 40th-order FD coefficients from the Taylor expansion method in the space domain.

Open in new tabDownload slide

Regarding the proposed method’s efficiency, we take the peak numerical error as a limitation to perform an experiment. In the homogeneous model with v = 2000 m s⁻¹, the grid system is 200 × 200, the spatial step is 10 m, the temporal step is 0.0002 and the propagation time is 0.4 s. The dominant 20 Hz frequency Ricker wavelet is located at (1000 × 1000 m). The corresponding result is illustrated in figure 8, showing that the proposed approach should adopt the 16th order FD coefficients and ADMM should use 20th order FD coefficients under the same peak numerical error. The 20th and 16th order FD coefficients involve 21 points and 17 points to compute the spatial derivative, respectively. Thus, the proposed approach attains a (21–17)/21 (19.05%) efficiency increase, compared to the ADMM.

Figure 8.

Trace comparison with the same peak numerical error for the ADMM and proposed approach. The reference solution is computed by the 40th-order FD coefficients from the Taylor expansion method in the space domain.

Open in new tabDownload slide

We also use the Marmousi model, where the spatial grid is 5 m, to avoid accumulating temporal errors; the temporal step is 0.0001 s and the grid system is 4078 × 2098. The injected source located at (10 195 m, 5 m) is the Ricker wavelet, whose dominant frequency is 20 Hz. The optimal coefficient is the FD 16th-order coefficients. Numerical results are shown in figures 9 and 10.

Figure 9.

The combined snapshots and wavefield errors for four methods in the homogeneous model. The first column refers to snapshots at 2, 3 abd 4 s, respectively. The second column is the residual error between the reference solution and snapshots calculated by four types of the optimal FD operator. The reference solution is computed by the 40th-order FD coefficients from the Taylor expansion method in the space domain.

Open in new tabDownload slide

Figure 10.

Traces from figures 7B, D and F, where 8A, 8B and 8C are the traces with x coordinate of 10 195 m from figures 7B, D and F, respectively. We highlight the specific peak error on the y-axis for all optimized methods.

Open in new tabDownload slide

Figure 9 infers that the snapshots of the proposed method have the smallest residual errors in the red box compared with those of the ADMM, the LS and the REA under the same error limitation. The single channel analysis in figure 10 demonstrates the proposed method’s advantage, as it affords the smallest peak error compared to the other methods. This experiment illustrates that even in the complex model simulation, the advantages of the proposed method concluded from the homogeneous modeling still exist.

5. Conclusion

The |${L}_1$| norm cost function obtains the smallest dispersion error within the bandwidth compared with the other norm cost function. In this paper, we use the dual ascent method to solve the Lagrange dual problem aim to obtain the optimal FD coefficients with a small dispersion error. In the proposed method, we exploit the |${L}_1$| norm model’s advantage to suppress of dispersion error and build on it the optimization function, construct the restricted condition by modified special M wavenumbers and then combine the optimization function and the restricted condition into the Lagrange dual problem. The optimization problem is solved by the dual ascent method. Under the same FD coefficient length and error limitation, the comparison of numerical experiments on the four methods verify that the optimal FD coefficients of the new method have the least dispersion error within the bandwidth, which is beneficial in the quick reduction of the accumulated errors and highly accurate simulation results for homogeneous and inhomogeneous models.

Data availability

The article’s data will be shared on reasonable request to the corresponding author.

Acknowledgements

This study is supported by the Major Scientific and Technological Projects of Shandong Energy Group (no. SNKJ2022A06-R23), the Marine S&T Fund of Shandong Province for Pilot National Laboratory for Marine Science and Technology (Qingdao) (no. LSKJ202203405), the Funds of Creative Research Groups of China (no. 41821002) and the Major Scientific and Technological Projects of CNPC (no. ZD2019-183-003).

Conflict of interest statement. None declared.

Appendix A. Proposed method based on absolute error

The absolute error of dispersion relation can be expressed as

Then we change equation (7) into

which is a new cost function based on the absolute error over the given range |$[ {0,{\beta }_{max}} ]$|⁠. After the same operations as equations (7)–(9), the linear system of the new cost function can be written as |$E\_in{t}_{absolute}( {{\bf c}} ) =\ \parallel\! {{\bf Ac}} + {{\bf b}}\! \parallel_1$|⁠, and A, b, c are as follows

|${\rm{w}}$|here |${{\rm{A}}}_{i,j}( {i = 1,2,3,\ \ldots \ N,j = 1,2,3,\ \ldots \ ,M} )$| represents the element in the coefficient matrix A, whose location is the ith row and the jth column, and |${{\rm{b}}}_i$| denotes the ith element in the vector b.

The new Lagrange dual problem with an absolute error cost function can be written as

The |${R}_{absolute}( {{\bf c}} ) = 0$| in the Lagrange dual problem is constructed by the modified special M wavenumbers of M−1 zero points. The process of altering the wavenumbers of M−1 zero points can be expressed as

The |$\beta \_zer{o}_i$| is the special M wavenumbers of the REA, where the cost function of the REA is also based on the absolute error of dispersion relation (equation (A.1)). The number of |$\beta \_zer{o}_i$| is M + 1, including the |$\beta \_zer{o}_1 = 0$|⁠, which is decided by properties of the equation (A.1) and not determined by the optimal FD coefficients. Hence, the |$\beta \_zer{o}_1 = 0$| can not be used in equation (A.6). For the new optimal problem based on absolute error, |$p = \frac{{\beta \_zer{o}_2}}{2}$|⁠, while for the optimization problem based on a relative error, |$p = \beta \_zer{o}_1$|⁠, because the special M wavenumbers cannot change too much, half of the first interval is appropriate. The excessive modification will result in the irregular oscillation of the dispersion curve, as illustrated in figure A1.

Figure A1.

Absolute error of the 12th-order optimal FD operators among the different values of p.

Open in new tabDownload slide

The modified wavenumbers restricted condition in the Lagrange dual problem can be written as

Solving process of the Lagrange dual problem based on absolute error is the same as 21–25. Within the same error limitation or effective bandwidth, the optimal FD coefficients of the new method with the absolute error cost function also have less of a low-wavenumber dispersion error than that of ADMM with the absolute error cost function, as shown in figure A2.

Figure A2.

Absolute error of the 12th-order optimal FD operators for the different optimal strategies.

Open in new tabDownload slide

References