有限差分数值模拟在油气藏的勘探开发过程中具有重要的作用.有限差分数值模拟可以帮助研究地震波在复杂介质中的传播规律和波场特征，对实际的勘探开发具有指导意义.但是有限差分法求解波动方程时由于利用离散网格点的差分代替时间与空间的偏导数，会出现数值频散现象，导致波场模拟精度的降低.为了消除波场模拟中的数值频散问题，许多专家学者在这方面做了大量的工作.Alford(Alford et al.，1974)在研究声波方程模拟精度时指出，对于时间2阶、空间4阶精度有限差分来说，当半功率点频率对应的一个波长内不少于5个网格点时能得到准确的结果.Dablain(Dablain，1986)的研究结果标明，对于时间4阶、空间10阶有限差分，Nyquist频率对应的一个波长内有3个网格点时可以得到有效的无频散的结果，通过提高阶数降低了单位波长内网格点的需求个数.Boris(Boris and Book，1973；Book et al.，1975)在研究流体运移时提出FCT通量校正法，并由Fei Tong(Fei and Larner，1995)将其引入到弹性波正演模拟中用于消除数值频散，基本原理为假设所有极值点都为数值频散引起，对所有网格点进行扩散校正处理，再对非局部极值点进行逆扩散校正.杜启振(杜启振等，2010)通过引入强约束条件和弱约束条件构造了不同的Lagrange函数，将求取条件极值得到的优化差分算子应用于VTI介质的数值模拟，有效的压制了数值频散.梁文权(梁文全等，2013)提出时间-空间域特定频率点满足频散关系的方法可以代替空间域泰勒展开法提高数值模拟的精度和效率.Liu Yang(Liu and Sen，2009b；Liu and Sen，2011)提出了在时间-空间域对波数进行泰勒展开，通过匹配系数确定有限差分系数，提高了时间-空间域有限差分算子精度.严红勇等(严红勇和刘洋，2013)采用一种基于时空域频散关系的有限差分方法来求解声波方程，频散分析和数值模拟表明时空域有限差分方法较传统有限差分方法模拟精度更高.李娜(李娜等，2014)结合最小二乘思想推导出新的频散改进差分系数，可以更有效地压制数值频散.在频率域正演方面，吴国忱(吴国忱和王华忠，2005)针对TI介质弹性波正演模拟指出可从对差分算子校正、提高差分方程的阶数和FCT方法三方面考虑压制数值频散.殷文(殷文等，2006)采用优化系数的25点优化差分算子比传统差分算子能更好的抑制数值频散.刘璐(刘璐等，2013)提出优化15点频率-空间域有限差分正演模拟，证实了此方法可以较好的压制频散.张衡(张衡等，2014)提出基于平均导数方法的声波方程频率域高阶正演，综合了常规网格和混合网格的算法，具有很好的精度和效率.

本文以双相介质弹性波方程的一阶应力-速度方程为例，阐述了交错网格高阶有限差分数值模拟方法，并通过最小二乘法将不同积分区间的频散曲线对比选取最佳积分区间.利用具有最佳积分区间的最小二乘法替代传统Taylor展开方法求取差分系数并进行数值模拟，结果表明，最小二乘法得到的优化差分系数既能保证模拟的准确性，又能比传统的Taylor展开法更好的压制数值频散，提高了有限差分数值模拟的精度.

各向同性双相介质二维弹性波动方程的一阶应力-速度方程(苏云，2010)形式可表示为

交错网格差分( Virieux，1984，1986；)的主要思想是定义主体网格和交错网格，并在上面定义差分方程的参数(包括质点速度、应力和双相介质弹性系数)，通过中心差分计算空间导数值，本文为十阶空间差分，即需要用到周围的十个点计算中间一个点的空间导数值.对于函数f(x)，其一阶空间导数写成2N阶交错网格差分格式(董良国等，2000a)可近似表达为

对于传统的Taylor级数展开法交错网格有限差分，差分系数为

本文采用最小二乘( Liu，2003)的方法对一阶空间导数微分与差分之间的误差积分，通过改变积分区间对比分析不同的积分区间所对应的频散曲线，选取一种最优化的积分区间，从而确定一系列最优化的差分系数.根据平面波理论，令

式(8)左端是一阶空间导数的微分保留项，右端是经过差分化简后的保留项.令

分别采用Taylor展开法和不同积分区间的最小二乘法求解空间一阶导数的差分系数并进行数值频散分析(Liu and Sen，2009a).图 1是空间差分阶数为10时，Taylor展开法以及积分上限分别为$\frac{\pi }{{16}}$，$\frac{{2\pi }}{{16}}$，…，$\frac{{8\pi }}{{16}}$的最小二乘法频散关系曲线.从图 1可以看出：(1)b=$\frac{\pi }{{16}}$，当β>0.4时，E(β)>1且逐渐增大，β＝1.4时数值频散最大，随后减小.(2)b= $\frac{{2\pi }}{{16}}$，当β>0.8时，E(β)<1且逐渐减小，数值频散逐渐增大.(3)b=$\frac{{3\pi }}{{16}}$，当β>0.85时，且E(β)<1逐渐减小，数值频散逐渐增大.(4)b= $\frac{{4\pi }}{{16}}$ 与b=$\frac{{5\pi }}{{16}}$，当β>0.9与β>1.0时，E(β)<1且逐渐减小，数值频散逐渐增大.(5)b=$\frac{{6\pi }}{{16}}$16，当β<0.25时，E(β)>1，当β>1.15时，E(β)<1且逐渐减小，数值频散增大.(6)b=$\frac{{7\pi }}{{16}}$，b=$\frac{{8\pi }}{{16}}$，以及Taylor展开，E(β)与精确值1偏离严重，且有较大的波动.综上所述，我们应该选择积分上限b=$\frac{{5\pi }}{{16}}$ 或者b=$\frac{{6\pi }}{{16}}$，具体应该依据双相介质模型中β的范围决定.

图 1

Fig. 1

图 1 数值频散分析 (a)总体频散曲线图；(b)局部频散曲线图. Fig. 1 The numerical dispersion analysis diagram (a)Total dispersion analysis diagram；(b)Partial dispersion analysis diagram.

稳定性条件是波动方程有限差分方法中重要的问题(侯安宁和何樵登，1995；孟凡顺等，2013).本文采用董良国(董良国等，2000b)提出的一阶弹性波方程交错网格高阶差分稳定性条件，对于二维各向同性介质，需要满足如下不等式：

经计算，当Δx=Δz时，二维各向同性介质稳定性条件为

为了验证基于双相介质交错网格有限差分最小二乘方法相对于传统Taylor展开方法的优越性，设计两个双相介质模型进行正演来模拟实际效果.

均匀双相介质模型，网格点数为500×500，网格间距为6 m×6 m，时间采样间隔为0.6 ms.震源为纵波震源，炮点坐标位于(1500 m，1500 m)，检波点纵坐标位置为1200 m.子波为雷克子波，主频是30 Hz.模型参数如表 1所示.分别采用积分上限b=5π/16的最小二乘法与Taylor展开方法的空间差分系数，时间精度为2阶，空间精度为10阶的交错网格差分格式进行正演模拟.图 2和图 3分别为两种方法480 ms时正演模拟的波场快照，图 4为波场快照中两种方法对频散压制效果的对比.图 5和图 6分别为两种方法相应的单炮地震记录，图 7为两种方法频散压制效果在地震记录中的对比.

表 1 (Table 1)

表 1 均匀双相介质弹性参数 Table 1 Elastic parameters of uniform two-phase medium 地层 A N Q R B ρ11 ρ12 ρ22 参数 12.72 6.84 0.953 0.331 3.0 2170 -83 191 注：A、N、Q、R的单位是109 kg·m-1·s-2，ρij的单位是kg·m-3，B的单位是N·s·m-4.

表 1 均匀双相介质弹性参数 Table 1 Elastic parameters of uniform two-phase medium

图 2

Fig. 2

图 2 均匀介质Taylor展开法弹性波t=480 ms时波场快照 (a)固相水平分量；(b)固相垂直分量；(c)流相水平分量；(d)流相垂直分量. Fig. 2 Elastic snapshots of Taylor series expansion method at t=480 ms in homogeneous medium (a)Solid and horizontal component；(b)Solid and vertical component； (c)Fluid and horizontal component；(d)Fluid and vertical component.

图 3

Fig. 3

图 3 均匀介质最小二乘法弹性波t=480 ms时波场快照 (a)固相水平分量；(b)固相垂直分量；(c)流相水平分量；(d)流相垂直分量. Fig. 3 Elastic snapshots of Least square method at t=480 ms in homogeneous medium (a)Solid and horizontal component；(b)Solid and vertical component； (c)Fluid and horizontal component；(d)Fluid and vertical component.

图 4

Fig. 4

图 4 均匀介质Taylor展开法与最小二乘法t=480 ms时的波场快照对比图 (a)Taylor展开法固相水平分量；(b)最小二乘法固相水平分量； (c)Taylor展开法流相垂直分量；(d)最小二乘法流相垂直分量. Fig. 4 Elastic snapshots of Taylor series expansion method and Least square method at t=480 ms in homogeneous medium (a)Solid and horizontal component of Taylor method；(b)Solid and horizontal component of Least square method； (c)Fluid and vertical component of Taylor method；(d)Fluid and vertical component of Least square method.

图 5

Fig. 5

图 5 均匀介质Taylor展开法弹性波正演模拟地震记录 (a)固相水平分量；(b)固相垂直分量；(c)流相水平分量；(d)流相垂直分量. Fig. 5 Elastic seismic record of Taylor series expansion method in homogeneous medium (a)Solid and horizontal component；(b)Solid and vertical component； (c)Fluid and horizontal component；(d)Fluid and vertical component.

图 6

Fig. 6

图 6 均匀介质最小二乘法弹性波正演模拟地震记录 (a)固相水平分量；(b)固相垂直分量；(c)流相水平分量；(d)流相垂直分量. Fig. 6 Elastic seismic record of Least square method in homogeneous medium (a)Solid and horizontal component；(b)Solid and vertical component； (c)Fluid and horizontal component；(d)Fluid and vertical component.

图 7

Fig. 7

图 7 均匀介质Taylor展开法与最小二乘法弹性波正演模拟地震记录对比图 (a)Taylor展开法固相水平分量；(b)最小二乘法固相水平分量； (c)Taylor展开法流相垂直分量；(d)最小二乘法流相垂直分量. Fig. 7 Elastic seismic record of Taylor series expansion method and Least square method in homogeneous medium (a)Solid and horizontal component of Taylor method；(b)Solid and horizontal component of Least square method； (c)Fluid and vertical component of Taylor method；(d)Fluid and vertical component of Least square method.

数值模拟结果表明：

1 )Taylor展开法正演模拟中弹性波固、流相水平分量和垂直分量均存在数值频散现象，模拟精度较低.并且慢纵波的频散现象相对快纵波较严重，这是由于慢纵波速度较低的结果.

2)最小二乘法正演模拟精度较高，快纵波与慢纵波的波场轮廓清晰，频散得到了很好的抑制.

设计一个凹陷模型，如图 8所示.网格点数500×500，网格间距为10 m×10 m，时间采样间隔为1 ms.震源为纵波震源，炮点坐标位于(2500 m，2500 m)，检波点纵坐标位置为2000 m.如图所示.子波为雷克子波，主频是30 Hz.模型参数如表 2所示.分别采用积分上限b=5π/16的最小二乘法与Taylor展开方法的空间差分系数，时间精度为2阶，空间精度为10阶的交错网格差分格式进行正演模拟.图 9和图 10分别为两种方法600 ms时正演模拟的波场快照，图 11为波场快照中两种方法的频散效果对比.图 12和图 13分别为相应的单炮地震记录，图 14为两种方法频散效果在地震记录中的对比.

图 8

Fig. 8

图 8 凹陷模型 Fig. 8 The sunken model

表 2 (Table 2)

表 2 凹陷模型弹性参数 Table 2 Elastic parameters of the sunken model 地层 A N Q R B ρ11 ρ12 ρ22 上层 12.72 6.84 0.953 0.331 3000 2170 -83 191 下层 31.65 12 0.7 0.79 5000 2430 -100 255 注：A、N、Q、R的单位是109 kg·m-1·s-2，ρij的单位是kg·m-3，B的单位是N·s·m-4.

表 2 凹陷模型弹性参数 Table 2 Elastic parameters of the sunken model

图 9

Fig. 9

图 9 凹陷模型Taylor展开法弹性波t=600 ms时波场快照 (a)固相水平分量；(b)固相垂直分量；(c)流相水平分量；(d)流相垂直分量. Fig. 9 Elastic snapshots of Taylor series expansion method at t=600 ms of sunken model (a)Solid and horizontal component；(b)Solid and vertical component； (c)Fluid and horizontal component；(d)Fluid and vertical component.

图 10

Fig. 10

图 10 凹陷模型最小二乘法弹性波t=600 ms时波场快照 (a)固相水平分量；(b)固相垂直分量；(c)流相水平分量；(d)流相垂直分量. ①凹陷底面反射快纵波；②凹陷底面反射慢纵波；③凹陷底面透射慢纵波；④凹陷底面透射快纵波； ⑤断面反射快纵波；⑥两侧界面透射快纵波；⑦两侧界面反射快纵波；⑧断面波经凹陷底面反射快纵波； ⑨断面波经凹陷底面透射快纵波；⑩震源直接激发快纵波； 震源直接激发慢纵波. Fig. 10 Elastic snapshots of Least square method at t=600 ms of sunken model (a)Solid and horizontal component；(b)Solid and vertical component； (c)Fluid and horizontal component；(d)Fluid and vertical component. ①Reflected fast P-wave of sunken underside；②Reflected slow P-wave of sunken underside；③Transmitted slow P-wave of sunken underside； ④Transmitted fast P-wave of sunken underside；⑤Reflected fast P-wave of section；⑥Transmitted fast P-wave of side face；⑦Reflected fast P-wave of side face；⑧Reflected fast P-wave of section wave via sunken underside；⑨Transmitted fast P-wave of section wave via sunken underside；⑩Reflected fast P-wave by source； Reflected slow P-wave by source.

图 11

Fig. 11

图 11 凹陷模型Taylor展开法与最小二乘法t=600 ms时的波场快照对比图 (a)Taylor展开法固相水平分量；(b)最小二乘法固相水平分量； (c)Taylor展开法流相垂直分量；(d)最小二乘法流相垂直分量. Fig. 11 Elastic snapshots of Taylor series expansion method and Least square method at t=600ms of sunken model (a)Solid and horizontal component of Taylor method；(b)Solid and horizontal component of Least square method； (c)Fluid and vertical component of Taylor method；(d)Fluid and vertical component of Least square method.

图 12

Fig. 12

图 12 凹陷模型Taylor展开法弹性波正演模拟地震记录 (a)固相水平分量；(b)固相垂直分量；(c)流相水平分量；(d)流相垂直分量. Fig. 12 Elastic seismic record of Taylor series expansion method of sunken model (a)Solid and horizontal component；(b)Solid and vertical component； (c)Fluid and horizontal component；(d)Fluid and vertical component.

图 13

Fig. 13

图 13 凹陷模型最小二乘法弹性波正演模拟地震记录 (a)固相水平分量；(b)固相垂直分量；(c)流相水平分量；(d)流相垂直分量. Fig. 13 Elastic seismic record of Least square method of sunken model (a)Solid and horizontal component；(b)Solid and vertical component； (c)Fluid and horizontal component；(d)Fluid and vertical component.

图 14

Fig. 14

图 14 凹陷模型Taylor展开法与最小二乘法弹性波正演模拟地震记录对比图 (a)Taylor展开法固相水平分量；(b)最小二乘法固相水平分量； (c)Taylor展开法流相垂直分量；(d)最小二乘法流相垂直分量. Fig. 14 Elastic seismic record of Taylor series expansion method and Least square method of sunken model (a)Solid and horizontal component of Taylor method；(b)Solid and horizontal component of Least square method； (c)Fluid and vertical component of Taylor method；(d)Fluid and vertical component of Least square method.

数值模拟结果表明：

1 )Taylor展开法正演模拟中凹陷模型的波场快照和地震记录存在十分严重的数值频散现象，对准确的分析波场特征存在着很大影响.

2)最小二乘法正演模拟波场快照与地震记录中快纵波的波形十分清晰，对分析复杂双相介质模型的波场特征具有重要作用.

3 )由于慢纵波速度较低，振幅衰减较快，t=600 ms时固相波场快照中震源激发的慢纵波衰减十分严重，波场几乎消失.但快纵波经过界面反射和透射转化的慢纵波清晰可见.固相介质中慢纵波的能量远小于流相介质中的慢纵波，且由于流体的耗散作用，慢纵波能量衰减很快，因此固相介质质点速度分量中很难观测到慢纵波.

4 )由于空间间隔取值较大，慢纵波速度很小，最小二乘法中慢纵波依然存在较严重的数值频散，如果需要可以通过进一步的降低网格间距、提高阶数等方法压制慢纵波的频散现象.

本文利用交错网格有限差分方法对均匀与凹陷双相介质模型进行了正演模拟.求取有限差分系数时，在最小二乘方法的基础上改变积分区间上限，针对不同的积分上限分别进行频散分析.选取压制频散效果最好的积分区间，求出最优化的有限差分系数并进行数值模拟.通过与传统Taylor展开方法的频散分析和数值模拟对比，这种优化差分系数的方法在不降低计算效率的情况下压制频散效果更好，精度更高.

通过基于优化差分系数的交错网格正演模拟方法与泰勒展开方法对比分析，得到如下几点结论：

1 )通过可变积分区间的最小二乘法对差分系数进行优化求取过程简单且容易实施.

2)优化差分系数法用于正演模拟效果明显，且符合传统方法的波场特征，在提高模拟精度的前提下并不会增加计算量.

3 )文中优化差分系数法可以通过改变网格间距、差分阶数等参数综合选择合适的模拟精度和模拟效率，满足自己的需要.

4)文中讲述的方法只是针对空间差分系数做了频散校正，对于时间差分以及更有效的压制时间与空间频散的方法有待研究.

致 谢 感谢国家高技术研究发展计划(973)(2013CB228604)多元多尺度地球物理信息映射关系及反演对本文的赞助.