因为声波方程有限差分方法计算效率高、所需内存相对较小、实现简单，而广泛应用于地震正演研究(Alford et al.，1974；Kelly et al.，1976；Basabe and Sen，2007；冯英杰等，2007；王珺等，2007；韩令贺等，2011；兰海强等，2011；Chu and Stoffa，2012a，2012b; Wang，2014)，同时也是地震逆时偏移成像技术得以快速发展的基础(李博等，2012;刘红伟等，2011，2012).

网格频散是有限差分中至关重要的问题，直接影响着有限差分法在波动方程数值解过程中的应用效果.网格频散是由对时间和空间偏导数网格化造成的，相速度变成了网格间距的函数(Dablain 1986)，这会导致地震波的数值相速度不等于地球介质的真实相速度，使得波场模拟的精度降低.一般来说，如果存在时间频散，则高频的相速度增大；如果存在空间频散，则高频的相速度减小(Dablain 1986).

为了缓解或者压制网格数值频散效应，一般采用两种措施：其一是采用低阶差分格式，要求时间步长和空间步长非常小，这会造成所需内存和计算量过大而无法应用于三维的情况；其二是采用高阶差分格式，一般情况下，时间方向采用二阶差分格式，空间方向采用高阶差分格式或伪谱法.高阶差分格式主要是利用等波纹准则或最小误差准则优化设计差分系数实现对空间偏导数的近似.伪谱法则是通过正、反傅里叶变换来实现空间偏导数的精确求解.两种方法相比，高阶空间差分格式具有精度适中、计算效率高的特征而得到广泛应用，特别是用于地震逆时偏移成像.伪谱法计算精度高，但因计算量过大，一般常用于地震正演模拟.此外，Dablain(Dablain 1986)和Chen(Chen 2007，2011)提出了时间方向导数的四阶有限差分格式，提高了波场模拟的精度.

在空间方向偏导数差分格式近似计算方面已开展了大量研究工作，积累了大量的研究成果.考虑到在地震波场数值模拟中，需在时间方向和空间方向同时网格化，两者共同决定着数值频散特征.因此仅依靠空间方向偏导数计算精度的提高，并不能压制时间频散，甚至加强了时间频散.鉴于上述原因，Finkelstein提出了在时间-空间域内设计与确定有限差分格式的新方法，以期提高正演模拟计算精度.其方法原理是空间域对波数进行泰勒展开、特定频率点满足相速度关系、群速度关系的组合，用以控制精度的阶数和频散(Finkelstein and Kastner，2006; Finkelstein and Kastner，2008).

偏微分方程的有限差分格式分为显式的和隐式的.显式有限差分格式中，某点导数的计算需要该点和它周围点的值.隐式有限差分格式中，某点导数的计算除了需要该点和它周围点的值之外，还需要它周围点的导数值(Liu and Sen，2009).使用隐式有限差分格式需要求解多对角阵，因而计算量较大.但是为了达到相同的精度，使用隐式有限差分算子比使用显式有限差分算子一般需要更少的网格点.Liu和Sen(2009)的研究指出，对于二阶导数来说，隐式差分格式(2N+2)阶有限差分算子的精度相当于(4N+2)阶显式有限差分算子的精度.但是传统的隐式有限差分有两个缺点：第一是在空间域进行泰勒展开求取有限差分系数，仅能在低频或者低波数较好的保持频散关系；第二是在空间域求取隐式有限差分系数，不能减小时间频散.因此，本文提出了使用线性化的方法在时间-空间域确定隐式有限差分算子，使得算子在更大的波数范围保持时间-空间域频散关系，从而提高正演模拟的精度和效率.

本文以一维声波方程为例分析差分格式的优化问题：

其中，p(x，t)是波场，v是纵波速度.

二阶导数的显式有限差分算子公式：

其中，p_mⁿ=p(x+mh，t+nτ)，h是空间步长，τ是时间步长.

函数p(x，t)二阶中心差分：

二阶导数的隐式有限差分算子公式为(Liu and Sen，2009)

其中，c_m(m=0，1，2，…，M)和b是隐格式有限差分系数.当b=0时，得到的是显式的二阶导数有限差分算子.

将公式(2)和(4)代入公式(1)，然后代入平面波p_mⁿ=p₀e^{i[k(x+mh)-ω(t+nτ)]}得到(Liu and Sen，2009)：

其中，r=vτ/h，α=kh/2，ω=kv.

公式(5)中，如果b=0，可以得到声波方程时间-空间域显式有限差分系数.如果b≠0，可以得到声波方程时间-空间域隐式有限差分系数.如果τ→0则

则时间-空间域隐式有限差分格式退化为空间域隐式有限差分格式.

根据震源最高频率，波速和网格间距确定需要优化的波数上限(梁文全等，2013)：

其中，f是最高震源频率.

在确定的波数范围内，选择M+1个均匀分布的波数点满足时间-空间域隐式差分格式有限差分算子频散关系(公式(5))，则得到下面的线性方程组：

解线性方程组，可以求得时间-空间域隐式有限差分算子的系数，其中

由公式(8)可见，时间-空间域隐式有限差分算子的系数与波速、所取波数范围、时间步长和空间步长有关.也就是说，时间-空间域中，对于不同的波速，隐格式有限差分的系数b是不同的.但是隐格式有限差分中，空间偏导数的求取需要解一个三对角阵，这个三对角阵主要有b构成，一般希望b值是固定的.同时我们知道，低波速需要的算子长度最长.因此根据最低波速确定b值，然后把b作为已知，求取系数c_m(m=0，1，2，…，M).空间域隐格式有限差分算子，空间域交错网格隐格式有限差分算子，时间-空间域隐格式交错网格有限差分算子的推导与上述相同，不再赘述.

使用下面的公式衡量时间-空间域频散(Liu and Sen，2009)：

δ的值越接近1，数值频散越小；否则数值频散越大.

由图 1可以看到，随着有限差分算子长度的增加，传统隐式有限差分算子和时间-空间域隐式有限差分算子的数值频散都减小.时间-空间域新方法的数值频散在整个kh范围内分布比较均匀，在更广的kh范围内保持较小误差，且可以根据震源频率、空间步长调整需要满足频散关系的波数上限使得kh比较小的情况下数值频散也很小.

图 1

Fig. 1

图 1 两种不同方法随M增加的数值频散图 (a) 频散对比； (b) a图的局部放大, 其中，v=2000 m·s-1时间步长h=10 m,空间步长τ=0.001 s. Fig. 1 Plot of dispersion curves for different M (a) the comparison of FD dispersion; (b) zoom of (a)

由图 2和3可以看出，随速度和时间步长的变化，传统隐式有限差分算子的数值频散变化比较大，而时间-空间域隐式有限差分算子频散误差变化比较小，因此新方法比传统的方法稳定.

图 2

Fig. 2

图 2 两种不同方法随M增加的数值频散图 (a) 频散对比； (b) a图的局部放大； 其中，v=2000 m·s-1时间步长h=10 m,空间步长τ=0.002 s Fig. 2 Plot of dispersion curves for different M when τ=0.002 s (a) the comparison of FD dispersion; (b) zoom of (a)

图 3

Fig. 3

图 3 两种不同方法随M增加的数值频散图 (a) 频散对比； (b) a图的局部放大；其中，v=4500 m·s-1时间步长h=10 m,空间步长τ=0.001 s Fig. 3 Plot of dispersion curves for different M when v=4500 m·s-1 (a) the comparison of FD dispersion; (b) zoom of (a)

首先使用均匀介质模型，模型大小为220×200.对于所有的隐式有限差分算子M=8，v=2000 m·s^-1，h=10 m，时间步长 dt=1 ms.震源放在中心位置. 震源子波定义为

其中，f₀=150 Hz，f₀/π是震源主频，t₀=4/f₀，c是常数.

传统隐式有限差分格式和时间-空间域隐式有限差分格式在350 ms时的波场快照如图 4所示.图 4的切片(z/dz=110)如图 5所示，我们以时间步长0.5 ms的伪谱法为解析解.由图 4可见，伪谱法(τ=1 ms)的时间频散严重，传统的隐式差分格式时间频散和空间频散都严重，时间-空间域隐式有限差分格式同时有效压制了时间频散和空间频散.

图 4

Fig. 4

图 4 传统隐式有限差分算子(a)和新时间-空间域有限差分算子(b)的波场快照 Fig. 4 Comparison of snapshots for the conventional implicit method and the new time-space domain implicit method

图 5

Fig. 5

图 5 图 4的切片图 Fig. 5 Slice of snapshots in Fig.4

图 6显示了Marmousi模型，其速度在1500 m·s^-1和4700 m·s^-1之间.震源子波如公式(11)所示，其中f₀=120 Hz.对于不同的有限差分算子h=10 m，dt=1 ms，M=8.

图 6

Fig. 6

图 6 Marmousi模型 Fig. 6 Marmousi velocity model

图 7a是空间域泰勒展开法隐式有限差分格式的地震记录，频散明显；图 7b是时间-空间域隐式有限差分算子的地震记录，时间和空间频散都得到压制.同时我们注意到，使用空间域隐式有限差分格式进行正演模拟的时候，在1000 ms时已经出现不稳定现象，数值模拟崩溃.而使用时间-空间域隐式有限差分格式则可以保持稳定.两种方法在x/dx=250处的地震记录如图 8所示，新方法明显同时减小了时间频散和空间频散.

图 7

Fig. 7

图 7 Marmousi模型地震记录 (a) 传统方法隐式有限差分算子; (b) 时间-空间域隐式有限差分算子. Fig. 7 Seismic records from the Marmousi model (a) The conventional implicit method; (b) The new time-space domain implicit method.

图 8

Fig. 8

图 8 Marmousi模型地震道对比 红线表示传统隐式方法计算的地震道，绿线表示时空域隐式方法计算的地震道. Fig. 8 Comparison of seismic seismograms from Marmousi model (the red line is the traditional implicit method, whereas the green line is the new time-space domain method).

本文提出了时间-空间域确定声波方程隐式有限差分算子的准则，即给定波数范围满足频散关系.方法本身无需在波数方向上进行泰勒展开，而是通过震源频率、波速和网格间距确定波数上限，这使得新的差分格式在更大范围符合波场传播规律.此外，与传统的隐式有限差分格式相比，新的隐式有限差分格式稳定性更好.

通过频散分析，以及二维的地震波数值模拟，验证了本文提出的新方法的有效性.因此，时间-空间域隐式有限差分算子方法可以代替空间域泰勒展开法隐式有限差分算子用于地震的正演、逆时偏移中，以提高地震波数值模拟的精度和效率.