2011, Vol. 54
地震波场的数值模拟技术是在已知地下介质结构和参数的情况下，通过理论计算来研究地震波在地下介质中的传播规律，得到人工合成地震记录的一种技术[1].其中有限差分法是一种常用模拟方法，它将波动方程中波场函数的空间导数和时间导数用相应空间和时间的差分来代替[2]，其精度和计算速度都相对较为理想[3].
常规有限差分方法使用单一尺度的网格步长对模型进行剖分，这种处理方法必然会导致局部采样过量影响计算效率或局部采样不足影响模拟精度[4， 5].为解决这个问题，Moczo(1989)[6]首次提出可变网格的思路，即不再使用固定网格步长进行数值模拟，网格步长可以根据模型横向速度梯度进行调整，对于横向变速剧烈的区域和小尺度地质体采用精细网格剖分，在介质属性变化平缓的区域使用大网格剖分，这种方法可以在满足模拟精度的前提下提高模拟效率.
可变网格有限差分模拟方法可以实现在空间各个方向同时改变采样步长，并相应改变局部时间步长.目前可见的数值试验中变网格有限差分算法的空间网格采样步长比可达到十倍[7]，对于非均匀缝洞介质来说还远远达不到精确模拟的要求.本文将变网格方法网格步长比推广至百倍以上，达到了模拟孔缝洞介质的精度要求，并改进了一般变网格方法网格剖分方式，使该方法的模拟精度和效率有了进一步的提高.
2 方法原理 2.1 基于PML边界条件的有限差分方程常规网格基于PML 边界条件的有限差分方程需要将PML 衰减因子加入声波方程，并用差分格式对波动方程的导数项进行逼近.
三维声波方程为:
|
(1) |
其中u为波场值，V为介质的纵波速度.将波场进行分解，有
|
(2) |
令
|
(3) |
将(3)式加入PML 衰减项并进行时域傅里叶变换:
|
(4) |
其中
|
(5) |
其中d(x)控制x方向上的衰减.整理并进行时域傅里叶逆变换可得PML 时域控制方程:
|
(6) |
一阶、二阶时间导数的中心差分格式分别为:
|
(7) |
其中dt为时间采样间隔，un1，un-11 ，un+11 分别表示当前时刻，前一时刻和下一时刻u1 的值.将(6)代入(5)式，整理可得:
|
(8) |
在x方向上，常规空间M阶精度差分格式为[9]:
|
(9) |
其中dx为x方向的采样步长，cm为差分系数.
y方向和z方向也有相似结论，基于PML边界的时间二阶、空间M阶精度声波方程常规差分格式为:
|
(10) |
其中d(x)，d(y)，d(z)分别为三个方向的衰减函数，dx，dy，dz为空间采样步长.
2.2 可变网格方法基于PML边界条件的有限差分方程可变网格有限差分方法的主要思想为:在需要使用小网格精细模拟的区域使用小网格对模型进行重采样计算，其他区域使用大网格计算.可变网格模型剖分的示意图如图 1.
|
图 1 空间可变网格剖分示意图 Fig. 1 Variable spatial grid meshing configuration |
常规差分格式是对称的，即参与计算的点相对求导点的距离两两相等.变网格算法使用不同尺度的空间网格剖分模型，如按照常规方式顺序取点计算变网格区域附近的导数势必产生误差.空间步长与差分算法的数值频散有关，大小网格步长的不同导致差分算法数值频散的不同，也使得大小网格的交界成为一个弱的人为反射界面[10， 11].
变网格交界处的过渡区各点需要重新推导差分格式以减小人为干扰.对于空间十阶精度的差分格式，需要修正差分格式的点有5个.过渡区的处理方法由图 2所示.星形符号标示的点为计算导数的点，图中同时标注了参与该点导数计算的点.修正的过渡区差分格式保持两边计算点距中心点对称的特征，以期减少人为干扰.
|
图 2 变空间步长过渡区域处理方式示意图 Fig. 2 Finite-difference schemes in the transition zone with different spatial steps |
过渡区的第n个点，在x方向上，空间M阶精度变网格过渡区的差分格式为[7]:
|
(11) |
则基于PML 边界的时间二阶空间M阶精度声波方程变网格过渡区的差分格式为:
|
(12) |
其中N为大小网格步长比.c′i为变网格差分格式的差分系数.
导数计算进入到小网格区以后使用小网格参数的常规差分格式.其中计算所需要的部分波场值是不存在的，此时需要通过对大网格点进行插值.由小网格区进入大网格区的过渡带也使用同样的处理方法.
空间采样间隔发生改变时，为了保证时域显式差分方法的稳定性，必须调整时间采样间隔，使得最小空间采样步长满足稳定性条件.若整个模型都采用统一的时间采样间隔，必然会造成大空间步长区域的时间过采样[4， 5].为了提高算法的模拟效率，需要在计算小网格区域时加密时间片循环，即在大网格区域采用大的时间步长，小网格区域采用小时间步长.
改变时间步长需要在大时间步长区设置一个过渡带，利用大时间步长求出过渡带区域中下一时刻的波场值，采用内插方法得到小时间步长计算所需的各时刻波场值.图 3为变时间步长的示意图.图中使用十字星标注的点为需要插值的点.由于时间差分算子精度小于空间差分算子精度，时间维度上的精度敏感性较空间低[10]，因此内插引入的误差不会造成明显的人为反射.
|
图 3 变时间步长过渡区处理方式示意图 Fig. 3 Finite-difference schemes in the transition zone with different temporal steps |
使用PML 条件需要在研究区域四周引入完全匹配层[8， 12].对于二维模型，大网格和小网格区使用公式(10)，变网格过渡区使用公式(12)进行递推模拟.令模型内部各点函数d(x)=0，d(z)=0，边界的四个角点d(x)≠0，d(z)≠0，速度V等于模型角点的速度，模型上下边界匹配层区域d(x)=0，d(z)≠0，速度V分别与模型上下边界速度相等，左右边界的匹配层d(x)≠0，d(z)=0，速度V分别与模型左右边界速度相等.调整吸收参数和匹配层厚度，波场传播到匹配层区域时不会产生边界反射，而是按照指数规律迅速衰减.
使用均匀介质模型测试PML 边界的效果，模型大小为1km×1km, 纵波速度为5000m/s, 使用40倍变网格算法，大网格空间步长为10m, 时间步长为1ms, 小网格空间和时间步长分别为0.25 m, 0.025ms.分别将变网格区置于模型中央和模型边缘测试PML 边界的吸收效果.震源为主频50 Hz的Ricker子波位于模型中心.测试结果(图 4)显示PML 边界有较好的吸收效果，在模型边界处没有明显的反射波.即使变网格区域位于模型边界处也可以保证边界的吸收效果.同时这个结果也测试了可变网格算法不会在变网格交界处产生人为干扰，保证了模拟精度.
|
图 4 PML 边界吸收效果.均匀模型140ms波场快照，图中标示了变网格区域 (a)变网格区域位于模型中央，硬边界；(b)变网格区域位于模型中央，PML边界；(c)变网格区域位于x方向边界，硬边界；(d)变网格区域位于x方向边界，PML边界；(e)变网格区域位于z方向边界，硬边界；(f)变网格区域位于z方向边界，PML边界. Fig. 4 140 ms snapshots of a homogenous media model show PML boundary condition absorbing effect.Grid step variation areas are marked on the figures (a) Grid step variation area locates in the centre of the model with hard boundary condition；(b) Grid step variation area locates in the centre of the model with PML boundary condition；(c) Grid step variation area locates on the edge of the model in x direction with hard boundary condition；(d) Grid step variation area locates on the edge of the model in x direction with PML boundary condition；(e) Grid step variation area locates on the edge of the model inzdirection with hard boundary condition；(f) Grid step variation area locateson the edge of the model inzdirection with PML boundary condition. |
改变空间网格步长所使用的过渡方法存在一定的局限性.对于空间十阶精度的差分算法，小网格区域至少要达到5个大网格长度才能满足过渡方法的要求.而小网格区两边的过渡都要满足要求，出于波场稳定和不引入人为反射干扰等因素的考虑，一般小网格区域的长度要达到10个大网格长度以上.对于较低倍数的变网格算法，这个条件相当容易满足.而变网格倍数达到数百倍至千倍以上时，5个大网格长度的区域对内存消耗也很大.微小尺度网格模拟的地质体一般规模较小，小网格区域过大会严重降低模拟效率.
数值试验表明，当大小网格步长比达到一千倍时，显式差分的稳定性条件很难得到满足.一般意义上时间采样步长近似正比于空间采样步长的稳定性关系不再成立，目前经数值试验证实的常规可变方法变步长的上限为800倍.
改变模型网格剖分方式可以改进该方法.改进的多级可变网格剖分示意图如图 5.
|
图 5 多级可变网格剖分示意图 Fig. 5 Multilevel grid variation meshing configuration |
改进的变网格剖分方法将模型中需要利用小网格计算的区域用内外两层小网格进行剖分，并通过两次变网格达到最终需要的模拟精度.内外两层变网格的倍数均可以为任意常数.上文提到的大小网格区域空间和时间步长的过渡区处理方法仍适用于这种网格剖分方式.嵌套的变网格剖分方式就使得我们可以通过低倍数的变网格算法来实现高倍的网格步长变化.此时内层小网格区域的长度只需是外层小网格步长的5 倍就可以满足变网格过渡的条件，从而使得高倍数变网格计算可以在更小的区域内进行.对于计算精度要求高，复杂地质体尺度小的情况，这种方法更为灵活，同时相对于常规变网格算法，这种方法也可以节约计算需要的内存.
数值试验表明，多级可变网格方法实现千倍以上的网格步长变化并不会引起差分方法稳定性的突变.只要满足低倍数变网格方法的稳定性条件，多级可变网格算法就是稳定的.这样，在内存条件允许的情况下可以考虑通过高倍变网格算法对毫米尺度的地质体进行精确模拟.
3 模型试算数值算例包括三个模型:溶洞模型、裂缝模型和生物礁模型.所有模型均使用基于声波方程的时间二阶、空间十阶精度差分算子进行数值模拟计算，差分方法均满足稳定性条件，数值频散也被控制在不影响模拟精度的范围内.
3.1 溶洞介质模型碳酸盐岩储集层一般为非均匀介质储集层[13， 14].而碳酸盐岩储集层多为缝洞型储集层.为了对缝洞型储层的地震响应进行精确模拟，网格剖分的尺度应当与地质模型中微小裂缝与孔洞的尺度相当，即网格的尺度需达到厘米数量级或毫米数量级.对整个地质模型均使用小网格剖分将消耗大量的内存，现有的计算条件很难满足计算要求.为了解决这个问题，可以采用变网格差分算法对缝洞介质进行模拟，在缝洞发育区使用小网格剖分，对小尺度地质体的地震响应进行精细刻画，在其他区域使用大网格，既保证计算精度，也提高了计算效率.
本文设计了一个小尺度溶洞模型和一个裂缝模型，证实使用高倍数变网格有限差分方法可以高效精确地模拟缝洞型介质.
图 6为溶洞模型示意图.模型大小为1km×1km.在均匀介质背景下设计了三个半径不同的溶洞.其半径从左至右依次为0.6cm, 2.5m, 12m.背景介质速度为5000m/s, 溶洞内部介质速度为2000m/s.使用多级可变网格剖分方法对模型进行离散.其中外层小网格范围为:横向380~640 m, 深度为400~600m.内层小网格范围为横向420~500m, 深度为460~540m.模型大网格的采样间隔为10 m, 外层小网格采样间隔为1m, 内层小网格采样间隔为0.1m时间采样间隔为1 ms, 0.1 ms, 0.01 ms.最大变网格倍数为100 倍.使用的震源子波为Ricker子波，主频为70Hz.使用空间采样步长为10 m, 时间采样步长为1ms的常规有限差分方法对同一模型使用同样的参数进行模拟，所得结果作为参照.
|
图 6 溶洞模型示意图 (a)大网格采样结果，只有最大溶洞在大网格上有样点.图中五角星标注炮点位置，方框标出了多级变网格区域；(b)外层小网格采样结果，较小溶洞均有少量样点，溶洞边缘不平滑；(c)内层小网格采样结果，较小溶洞得到了较为精确的刻画. Fig. 6 Geometryofthecavitymodel (a) Sampling result with the largest spatial grid step.Only the biggest cave could be seen on the result.The shot point and the grid step variation domains are marked in the figure；(b) Sampling result with smaller grid step of the outer grid variation area.Two smaller caves could be sampled with ladder edge；(c) Sampling result with the smallest grid step of the inner grid variation area.All the caves could be sampled accurately. |
图 7为两种方法得到的单炮记录，为突出溶洞的响应切除了炮记录中的直达波.通过两种方法结果的对比可以看到，常规有限差分方法得到的只是最大溶洞的响应，这是由于常规差分算法的空间采样步长为10m, 两个较小溶洞的直径均小于采样步长，大网格离散后的模型只包含最大溶洞的样点.溶洞内部介质的速度较低，常规差分算法得到的炮记录上反映的不仅是溶洞的散射波，同时也伴随着很强的数值频散.变网格算法在溶洞区使用小网格进行采样，对于较大溶洞来说，增加采样点可以有效地压制由于空间采样步长过大造成的差分方法的数值频散，从而得到溶洞的正确的响应，溶洞散射波在炮记录上有清晰的体现，溶洞底部起到了汇聚能量的作用，可以看到能量明显增强的散射波，同时溶洞造成的多次波在炮记录上也有清晰的体现.两个较小的溶洞由于其尺度远小于地震子波的波长，在炮记录上的响应主要表现散射波.这两种方法的对比体现出了可变网格有限差分方法的精度要优于常规有限差分方法的精度.
|
图 7 溶洞模型单炮记录(已去直达波) (a)多级可变网格算法得到的单炮记录；(b)常规有限差分算法得到的单炮记录. Fig. 7 Shot records (without direct waves)of the cavity model (a) Shot record calculated by multilevel grid step variation method；(b)Shot record calculated by standard difference scheme. |
图 8为裂缝模型示意图.模型的大小为100m×100m, 在深度为50m 处有一个层界面，上层速度为4500m/s, 下层速度为6000 m/s, 界面中央有一个高度为3m 的小型突起，在突起界面处存在多条裂缝，裂缝宽度为5 mm, 密度为9 条/m, 方位角为90°，裂缝内填充介质的纵波速度为2000 m/s.对这个模型使用常规高倍可变网格算法，空间大网格采样步长为1m, 小网格步长为0.25mm, 较大时间步长为0.08ms, 较小时间步长为0.0002 ms.变网格区域的横向和深度范围都在46~54 m 之间.同时使用空间采样步长为1m, 时间采样步长为0.1ms的常规有限差分算法计算同样的模型，得到的结果与变网格算法的结果对比.两种方法使用的震源子波均为Ricker子波，子波主频为40Hz.
|
图 8 裂缝模型示意图 (a)大网格采样结果，不包含裂缝信息.图中五角星标注炮点位置，方框标注变网格区域；(b)为小网格采样结果，凸界面形状平滑，裂缝信息较完整. Fig. 8 Geometry of the fracture model (a) Sampling result with larger spatial grid step.The result does not contain any information of the fractures.The shot point and the grid step variation domain are marked in the figure；(b) Sampling result with smaller spatial grid step.The interface is described smoothly and the result contains complete fracture information. |
两种方法对裂缝模型的正演结果如图 9.为了使结果更清晰，炮记录均切除直达波，并增加10%的振幅增益.常规网格离散后的模型中没有反映裂缝的样点，因此在常规网格模拟方法得到的炮记录上只有水平界面的反射波与小型凸界面的反射，凸界面的反射波顶部明显上凸，其下层界面的反射波能量被凸界面反射波影响稍有减弱.变网格模拟方法在裂缝区使用了毫米数量级的采样间隔，确保了每个裂缝在离散后的模型上均有样点.在变网格算法得到的单炮记录上可以看出，由于裂缝的反射极性与地层的反射极性相反，裂缝的存在使得凸界面的反射波能量有所减弱.分布规则的裂缝发育带在炮记录上的响应与弱反射层面类似，其响应并非为单一的反射波，而是在一个区域内部存在大量散射波，分布于小型凸界面的散射波之下.
|
图 9 裂缝模型单炮记录(已去直达波) (a)可变网格算法得到的单炮记录；(b)常规有限差分算法得到的单炮记录. Fig. 9 Shot records( without direct waves)of the fracture model (a) Shot record calculated by grid step variation method；(b) Shot record calculated by standard difference scheme. |
裂缝模型证实变网格算法可以达到毫米级别的模拟精度.运行这个程序所需的内存为1G 左右，若为了精确模拟裂缝而对整个模型都使用小网格采样，则需要的内存为2T 左右，目前的计算机很难满足计算的要求，耗费的机时也相当可观.变网格算法在很大程度上提高了模拟效率，节约了内存和机时.
3.3 生物礁模型这里用一个我国西部某地区的生物礁模型来验证变网格方法的优越性.生物礁模型示意图如图 10.模型大小为5km×2km, 从上到下各层岩性和其纵波速度均在表 1中列出.生物礁被断层穿过，边缘不平滑，且生物礁内部分布着多个半径小于5m的溶洞，同时使用常规有限差分算法和40倍可变网格算法计算该模型，其中常规有限差分方法空间采样间隔为10m, 时间采样间隔为0.5 ms, 变网格算法空间大网格采样间隔为10m, 小网格采样间隔为0.25m, 时间采样间隔分别为0.5ms和0.0125ms, 变网格范围主要针对生物礁及其中的小型溶洞.横向范围为3.2~4.2km, 深度范围为500~900 m.两种方法使用的震源子波均为主频为40Hz的Ricker子波.
|
图 10 生物礁模型示意图( a)大网格采样结果，无生物礁内溶洞的信息，图中五角星符号标注炮点位置，方框标注变网格区域；(b)变网格区小网格采样结果，小尺度溶洞信息较为完整. Fig. 10 Geometry of the reef model (a)Sampling result with larger spatial grid step.The result does not contain any information of the caves distributed in the reef.The shot point and the grid step variation domain are marked in the figure；(b)Sampling result with smaller spatial grid step.Fine scale caves are described more accuracy. |
|
表 1 生物礁模型参数 Table 1 Parameters of the reef model |
图 11显示了两种不同方法得到的模拟结果.炮记录均增加10%的振幅增益.生物礁底部和边缘有断层分布，造成了生物礁体反射波伴随着很多断点的绕射波，常规有限差分方法和变网格方法均可以得到各层界面的响应和生物礁断点的复杂绕射波.生物礁内部有较多小型溶洞发育，常规有限差分方法使用大网格对模型进行采样，没有反映小型溶洞的样点，因此在常规有限差分方法得到的炮记录上没有生物礁内部溶洞的响应.变网格方法使用小网格对生物礁内部进行精细采样，较为精确地刻画了溶洞的特征，在变网格方法的炮记录上可以看到，在生物礁顶部反射波之下存在着大量的杂乱绕射波，其能量小于层面的反射.两种方法得到的炮记录的差波场可以清晰地看到多个小型溶洞的绕射波.溶洞相互影响使得其波场更为复杂.生物礁内部溶洞波场和其边缘断点绕射波使得波场中出现了一个地震模糊带，影响生物礁之下的层界面的反射波的振幅能量.
|
图 11 生物礁模型单炮记录 (a)可变网格算法得到的单炮记录；(b)常规大网格得到的单炮记录；(c)生物礁内部小型溶洞的响应(差波场). Fig. 11 Shot records of the reef model (a) Shot record calculated by grid step variation method；(b) Shot record calculated by standard difference scheme；(c) The subtraction of two wavefields. |
本文将一般变空间时间步长有限差分模拟方法推广到大小网格步长比为数百倍的高倍数变网格有限差分算法，并改进了一般变网格算法的网格剖分方式，使得变网格方法更灵活高效.经过数值试验验证该方法不会引入过多人为误差，并且高倍数的变网格算法可以精确有效地模拟复杂的碳酸盐岩孔缝洞型介质.
对于千倍以上的变网格算法，当小网格的尺度到达毫米的数量级时，所需的内存和计算机时的消耗均很大，程序的并行化是其发展方向之一.
| [1] | 冯英杰, 杨长春, 吴萍. 地震波有限差分模拟综述. 地球物理学进展 , 2007, 22(2): 487–491. Feng Y J, Yang C C, Wu P. The review of the finite-difference elastic wave motion modeling. Progress in Geophysics (in Chinese) (in Chinese) , 2007, 22(2): 487-491. |
| [2] | Dablain M A. The application of high-order differencing to scalar wave equation. Geophysics , 1986, 51: 54-66. DOI:10.1190/1.1442040 |
| [3] | 张永刚. 地震波场数值模拟方法. 石油物探 , 2003, 42(2): 143–148. Zhang Y G. On numerical simulations of seismic wavefield. Geophysical Prospecting for Petroleum (in Chinese) (in Chinese) , 2003, 42(2): 143-148. |
| [4] | Tessmer E. Seismic finite-difference modeling with spatially varying time step. Geophysics , 2000, 65: 1290-1293. DOI:10.1190/1.1444820 |
| [5] | Falk J, Tessmer E, Gajewski D. Efficient finite-difference modeling of seismic finite-difference modeling of seismic waves using locally adjustable time step sizes. Geophys| Prosp| , 1998, 46: 1290-1293. |
| [6] | Moczo P. Finite-difference technique for SH waves in 2D media, using irregular girds: application to the seismic response problem. Geophys| J| Int , 1989, 99: 321-329. |
| [7] | 黄超, 董良国. 可变网格与局部时间步长的高阶差分地震波数值模拟. 地球物理学报 , 2009, 52(1): 176–186. Huang C, Dong L G. High-order finite-difference method in seismic wave simulation with variable grids and local time-steps. Chinese J| Geophys| (in Chinese) (in Chinese) , 2009, 52(1): 176-186. |
| [8] | 王守东. 声波方程完全匹配层吸收边界. 石油地球物理勘探 , 2003, 38(1): 31–34. Wang S D. Absorbing boundary condition for acoustic wave equation by perfectly matched layer. OGP (in Chinese) (in Chinese) , 2003, 38(1): 31-34. |
| [9] | 刘洋, 李承楚, 牟永光. 任意偶数阶精度有限差分法数值模拟. 石油地球物理勘探 , 1998, 33(1): 1–10. Liu Y, Li C C, Mou Y G. Finite-difference numerical modeling of any even-order accuracy. OGP (in Chinese) , 1998, 33(1): 1-10. |
| [10] | 李振春, 张慧, 张华. 一阶弹性波方程的变网格高阶有限差分数值模拟. 石油地球物理勘探 , 2008, 43(6): 711–716. Li Z C, Zhang H, Zhang H. Variable-grid high-order finite-difference numeric simulation of first-order elastic wave equation. OGP (in Chinese) (in Chinese) , 2008, 43(6): 711-716. |
| [11] | 吴国忱, 王华忠. 波场模拟中的数值频散分析与校正策略. 地球物理学进展 , 2005, 20(1): 58–65. Wu G C, Wang H Z. Analysis of numerical dispersion in wave-field simulation. Progress in Geophysics (in Chinese) (in Chinese) , 2005, 20(1): 58-65. |
| [12] | Berenger J P. A perfectly matched layer for the absorption of electromagnatics waves. J|Comput| Phys| , 1994, 114: 185-200. |
| [13] | 刘远志. 碳酸盐岩地震相分析与数值模拟. 成都:成都理工大学, 2009. Liu Y Z. Seismic facies analysis to carbonate and numerical simulation . Chengdu: Chengdu University of Technology,2009 |
| [14] | 吴俊峰, 姚姚, 撒利明. 碳酸盐岩特殊孔洞型构造地震响应特征分析. 石油地球物理勘探 , 2007, 42(2): 180–185. Wu J F, Yao Y, Sa L M. Analysis on seismic response of special cavernous structure of carbonate. OGP (in Chinese) (in Chinese) , 2007, 42(2): 180-185. |