The first-order velocity-stress acoustic wave-equation can be described as follows:

$$\frac{{\partial p}}{{\partial t}}= - \rho {v^2}\left( {\frac{{\partial {V_z}}}{{\partial z}}+\frac{{\partial {V_x}}}{{\partial x}}} \right),$$

1

$$\frac{{\partial {V_x}}}{{\partial t}}= - \frac{1}{\rho }\frac{{\partial p}}{{\partial x}},$$

2

$$\frac{{\partial {V_z}}}{{\partial t}}= - \frac{1}{\rho }\frac{{\partial p}}{{\partial z}}.$$

3

where *p* is the pressure, *ρ* is the density and *v* is the wave propagation speed, *V**x* and *V**z* are the particle velocities. A constant density is assumed for the equations in this paper.

The acoustic wave-equation (1-3) can be discretized as

$$p_{{i,j}}^{k}=p_{{i,j}}^{{k - 1}} - {v^2}\tau \left( {{\Delta _z}{V_z}_{{(i{\text{+}}1/2,j)}}^{{k{\text{+}}1/2}}+{\Delta _x}{V_x}_{{(i,j{\text{+}}1/2)}}^{{k{\text{+}}1/2}}} \right),$$

4

$${V_x}_{{(i,j+1/2)}}^{{k+1/2}}={V_x}_{{(i,j+1/2)}}^{{k - 1/2}} - \tau \left( {{\Delta _x}p_{{i,j}}^{k}} \right),$$

5

$${V_z}_{{(i+1/2,j)}}^{{k+1/2}}={V_z}_{{(i+1/2,j)}}^{{k - 1/2}} - \tau \left( {{\Delta _z}p_{{i,j}}^{k}} \right),$$

6

whereis the time step, \(\varsigma _{{m,j}}^{n}=\varsigma (z+mh,x+jh,t+n\tau ),\varsigma =[{V_z},{V_x},p]\). The symbol \(\varDelta\) represents the implicit discrete form of the spatial derivative, for example

$$q_{{i,j+1/2}}^{k}={\Delta _x}p_{{i,j+1/2}}^{k}=\frac{{\frac{1}{h}\sum\limits_{{m=1}}^{M} {{c_m}(p_{{i,j+m}}^{k} - p_{{i,j - m+1}}^{k})+\frac{1}{h}{c_{M+1}}(p_{{i+1,j+1}}^{k} - p_{{i+1,j}}^{k}+p_{{i - 1,j+1}}^{k} - p_{{i - 1,j}}^{k})} }}{{1+b{h^2}\frac{{{\delta ^2}}}{{\delta {x^2}}}}},$$

7

where

$$\frac{{{\delta ^2}q}}{{\delta {x^2}}}=\frac{{q(x+h)+q(x - h) - 2q(x)}}{{{h^2}}},$$

8

.

in which, *M* is the length of the SGFD operators, *c**m* and *b* are the SGFD coefficients and *h* is the grid size. *M*=7 and \({c_{M+1}} \ne {\text{0}}\) are assumed in this paper in order to use a larger time and space grid interval.

Then we can implicitly get an approximation to the spatial derivative \({{q \approx \partial p} \mathord{\left/ {\vphantom {{q \approx \partial p} {\partial x}}} \right. \kern-0pt} {\partial x}}\)as

$$bq_{{i,j+3/2}}^{k}+(1 - 2b)q_{{i,j+1/2}}^{k}+bq_{{i,j - 1/2}}^{k}=\frac{1}{h}\left( {\sum\limits_{{m=1}}^{M} {{c_m}(p_{{i,j+m}}^{k} - p_{{i,j - m+1}}^{k})+{c_{M+1}}(p_{{i+1,j+1}}^{k} - p_{{i{\text{+}}1,j}}^{k}+p_{{i - 1,j+1}}^{k} - p_{{i - 1,j}}^{k})} } \right).$$

9

Hereafter, we denote Eqs. (4), (5) and (6) as the ISGFD scheme. One can use a larger time step with the coefficient *c**M*+1. At the same time, one can use a short operator length with the coefficient *b*11, 12. We propose to use the simplest explicit second-order SGFD operator for the spatial derivatives in Eqs. (2) and (3). This approach differs from the previous SGFD methods that use the same ISGFD operator for all of the first-order spatial derivatives in Eqs. (1), (2) and (3). The proposed HEI-SGFD scheme is given by:

$$p_{{i,j}}^{k}=p_{{i,j}}^{{k - 1}} - {v^2}\tau \left( {{\Delta _z}{V_z}_{{(i{\text{+}}1/2,j)}}^{{k{\text{+}}1/2}}+{\Delta _x}{V_x}_{{(i,j{\text{+}}1/2)}}^{{k{\text{+}}1/2}}} \right),$$

10

$${V_x}_{{(i,j+1/2)}}^{{k+1/2}}={V_x}_{{(i,j+1/2)}}^{{k-1/2}}-\frac{\tau }{h}(p_{{i,j+1}}^{k}-p_{{i,j}}^{k}),$$

11

$${V_z}_{{(i+1/2,j)}}^{{k+1/2}}={V_z}_{{(i+1/2,j)}}^{{k-1/2}}-\frac{\tau }{h}(p_{{i+1,j}}^{k}-p_{{i,j}}^{k}).$$

12

The SGFD scheme in Eq. (10) is implicit while the SGFD scheme in Eqs. (11) and (12) are explicit; therefore, we refer the proposed SGFD scheme as the HEI-SGFD scheme. The HEI-SGFD scheme is an ISGFD scheme although some of spatial derivatives are approximated with explicit SGFD scheme. It seems that Eqs. (11) and (12) are not accurate. However, this is not the case. The SGFD coefficient in Eq. (10) is determined by considering Eqs. (11) and (12). It is easily observed that the proposed HEI-SGFD scheme in Eqs. (10)- (12) needs less floating-point computations compared with the ISGFD scheme given in Eqs. (4)-(6). The computational resources needed by Eqs. (11) and (12) are smaller than the computational resources needed by Eqs. (5) and (6).

Substitute Eqs. (11) and (12) into Eq. (10) and we get:

In theory, the first-order HEI-SGFD scheme in Eqs. (10)- (12) is equivalent to the second-order implicit FD scheme with ordinary grid format as shown in Eq. (14).

The dispersion relation of the proposed HEI-SGFD scheme is

The SGFD coefficient in Eq. (15) should be carefully determined with the optimization methods. Similar to the method proposed by Wang et al.20, we get the objective function from Eq. 15 by minimizing the error between the true velocity and the phase velocity:

$$\Phi (c)=\sum\limits_{{k={\text{3e-12}}}}^{K} {\sum\limits_{{\theta =0}}^{{\pi /4}} {{{\left( {\frac{{{v_{fd}}}}{v} - 1} \right)}^2}} } =\sum\limits_{{k={\text{ 3e-12}}}}^{K} {\sum\limits_{{\theta =0}}^{{\pi /4}} {{{\left( {\frac{{a{\text{cos}}\left( {1+\frac{{{r^2}g}}{{2\left[ {1+2b\left( {\cos \left( {{k_z}h} \right) - 1} \right)} \right]\left[ {1+2b\left( {\cos \left( {{k_x}h} \right) - 1} \right)} \right]}}} \right)}}{{kv\tau }} - 1} \right)}^2}} } ,$$

16

where

$$\begin{gathered} {\text{ }}g{\text{= }}\left[ { - 2\sin \left( {0.5{k_x}h} \right)\sum\limits_{{m=1}}^{M} {2{c_m}} \sin \left( {(m - 0.5){k_x}h} \right)+4{c_{M+1}}\cos ({k_z}h)(\cos ({k_x}h) - 1)} \right]\left[ {{\text{1+2}}b{\text{(cos(}}{k_z}h{\text{)}} - {\text{1)}}} \right] \hfill \\ {\text{ }}+\left[ { - 2\sin \left( {0.5{k_z}h} \right)\sum\limits_{{m=1}}^{M} {2{c_m}} \sin \left( {(m - 0.5){k_z}h} \right)+4{c_{M+1}}\cos ({k_x}h)(\cos ({k_z}h) - 1)} \right]\left[ {{\text{1+2}}b{\text{(cos(}}{k_x}h{\text{)}} - {\text{1)}}} \right]. \hfill \\ \end{gathered}$$

17

The only unknowns in Eq. (16) are \({c_m}\) (\(m=1, \cdots ,M+1\)) and *b*. The wave number *k* in Eq. (16) should start from zero. However, there is *k* in the denominator, which causes instability. Therefore, we let *k* start from a very small number, e.g., 3e-12. We assume that the parameters such as the wave propagation speed, the time step and the spatial grid intervals are already given (then *r* = \(v\tau\)/*h* will be known). The other two parameters are *k* and *θ*. The propagation angle *θ* is from 0 to *π*/4 (0, *π*/16, 2*π*/16 ... *π*/4). From Eq. (16), we found that in the frequency-wave number domain the SGFD coefficients are related to the Courant ratio *r*. When one considers a wave-equation simulation with both fixed time and spatial steps grid intervals, the SGFD coefficient is different for different velocities and the stencil forms a big 3D matrix for a 2D complex velocity model. In order to get the hybrid explicit-implicit SGFD coefficient we apply the MATLAB function *lsqnonlin* to solve the nonlinear least-squares problem in Eq. (16).

The proposed HEI-SGFD scheme can be easily extended to 3D:

$$p_{{l,m,n}}^{k}=p_{{l,m,n}}^{{k - 1}} - {v^2}\tau \left( {{\Delta _z}{V_z}_{{(l,m,n)}}^{{k - 1/2}}+{\Delta _x}{V_x}_{{(l,m,n)}}^{{k - 1/2}}+{\Delta _y}{V_y}_{{(l,m,n)}}^{{k - 1/2}}} \right),$$

18

$${V_x}_{{(l,m+1/2,n)}}^{{k+1/2}}={V_x}_{{(l,m+1/2,n)}}^{{k-1/2}}-\frac{\tau }{h}(p_{{l,m+1,n}}^{k}-p_{{l,m,n}}^{k}),$$

19

$${V_y}_{{(l,m.n+1/2)}}^{{k+1/2}}={V_y}_{{(l,m.n+1/2)}}^{{k-1/2}}-\frac{\tau }{h}(p_{{l,m,n+1}}^{k}-p_{{l,m,n}}^{k}),$$

20

$${V_z}_{{(l+1/2,m,n)}}^{{k+1/2}}={V_z}_{{(l+1/2,m,n)}}^{{k-1/2}}-\frac{\tau }{h}(p_{{l+1,m,n}}^{k}-p_{{l+1,m,n}}^{k}).$$

21

whereis the time step, \(\varsigma _{{l,m,n}}^{k}=\varsigma (z+lh,x+mh,y+nh,t+k\tau ),\varsigma =[{V_z},{V_x},{V_y},p].\)The SGFD coefficient in Eq. (18) of the hybrid explicit-implicit SGFD scheme for 3D can be determined similarly.