变网格有限差分d性波方程数值模拟

朱生旺1，2　曲寿利1　魏修成1　刘春园3

（1.中国石化石油勘探开发研究院，北京100038；2.中国地质大学（北京），北京100083；3.中国石油大学（北京），北京102249）

摘要 研究复杂介质中地震波传播规律及地震响应特征需要在细小网格剖分下进行d性波方程数值模拟计算，而小网格下的数值模拟计算将带来巨大的计算量问题，采用变网格计算是减少计算量的有效途径。本文给出一种变网格差分计算的实现方法，在局部复杂介质区域采用细网格计算，其余区域采用粗网格计算，在两种网格的过渡区通过改变差分算子和波场插值实现波传播的过渡衔接。该实现方法的特点是：①在所有计算节点均采用交错网格计算方式；②估算空间一阶导数时，保持差分计算节点的对称性；③可选择任意阶计算精度。该方法的计算结果表明，粗细网格的过渡不会对地震波传播模拟带来影响，从而达到了既减少计算量又保证计算精度的目的。

关键词 d性波方程 数值模拟 变网格差分 频散

Elastic Modeling on A Staggered Grid with Varying Spacing

ZHU Sheng-wang1，2，QU Shou-li1，WEI Xiu-cheng1，LIU Chun-yuan3

（1.Exploration ＆ Production Research lnstitute，SlNOPEC，Beijing100083；2.China University of Geoscience，Beijing100083；3.China University of Petroleum，Beijing102249）

Abstract It is necessary to use small grid in elastic wave equation modeling for researching the law of seismic wave propagation and feature of seismic response in complex medium.It is an effective scheme to calculate on a staggered grid with varying spacing which can reduce the enormous calculation amount caused by using the small grid calculation means.This article proposes the finite difference method in which small grid is used for local complex medium，and big grid is used for others.In transition region between two kinds of grid，transient conjugation of wave propagation is implemented by changing the difference operator and applying wave field interpolation.The trait of this way includes：（1）Keeping staggered-grid calculation mode in all operation node；（2）Maintaining The symmetry of operation node in estimating spatial first derivate；（3）Any precision order available in finite difference calculation.The result of the way indicates that the simulation of seismic wave propagation will not be affected by using two kinds of grid，and as a result，not only the calculation amount decreases，but also the high calculation accuracy remains.

Key words elastic wave equation modeling staggered-grid dispersion

在复杂油气藏的勘探开发中，地震数值模拟（正演）发挥着越来越重要的作用。利用地震方法进行储层预测，首先要研究储层的地震响应特征，找出与储层最为相关的地震属性，进而优选储层预测的技术方法，地震数值模拟是研究储层地震响应特征不可缺少的手段。

对复杂储层或地质体进行地震数值模拟需在足够小的空间网格上进行计算，以保证地质模型的局部变化细节能够准确地反映到数值模拟结果中来，即确保数值模拟结果能够较为精确地反映介质的小尺度非均质性引起的波场变化。在正演计算中，为保证计算收敛，空间和时间方向的采样率必须满足收敛条件，由d性波方程有限差分正演的收敛条件［1，2］知，减小空间计算网格尺度的同时，时间方向的递推计算步长也必须大体上做相应比例的减小，即空间计算网格点数增加 n倍，意味着计算量将增加 n2倍。因此，若对整个模型采用细小网格计算所面临的问题是计算量太大，即使是二维模型的正演计算，譬如使用1m×1m 或更小尺度的计算网格，当采用能够反映实际地震反射深度的地质模型时，模拟多次覆盖采集，其计算量亦十分惊人。实际中的数值模拟一般是针对某一特定的地质目标（或储层），且地质目标仅占据整个模型很小的一部分，对计算网格仅要求在该地质目标处细化，于是，采用变网格计算，即在目标地质体区域采用较小的计算网格，而在该区域以外采用较大的计算网格就可以大大地减少计算网格节点的数量，从而可在很大程度上降低计算量。

为减少计算量，提高数值模拟的计算效率，Jastram C和Tessmer E提出了单方向上变网格计算方法［3］，文献［2］亦强调采用变网格计算技术，但没有给出具体实现方法。采用变网格进行差分计算的关键是解决好网格尺度变化处的波场过渡衔接问题，保证在过渡区不会因网格尺度变化产生明显的计算噪声。本文基于二维d性波方程数值模拟的交错网格差分算法，通过引入不等间距节点情况下一阶导数的高阶精度差分计算式，结合波场的高阶精度插值和粗细网格过渡区的计算节点选取技巧，给出d性波方程变网格正演计算的实现方法，应用结果表明，该方法能够得到满意的计算精度，在过渡区不会产生明显的因网格尺度变化带来的附加计算噪声。

1 利用对称任意间距节点估算一阶导数的高阶精度差分近似式

在沿时间方向的递推计算过程中，提高波场函数空间一阶导数的估计精度是d性波方程数值模拟研究的关键内容之一。对于固定差分网格，通过应用等距节点高阶精度的差分算子可实现空间导数的高精度估算［1，2］。对于局部采用细网格的非固定差分网格情况，在粗、细网格区域内部可分别应用上述等距节点的高阶精度的差分算子，而在粗、细网格过渡区，本文采用下面导出的具有任意节点间距的差分算子估算一阶导数值。

设在x方向上以点x0为中心对称分布2N个网格节点，其x坐标分别为x0-qNΔx/2，…，x0-q1Δx/2，x0+q1Δx/2，…，x0+qNΔx/2，如图1 所示，这里Δx为节点间的最小间距。由2N个节点处的函数值f（x0-qNΔx/2），…，f（x0-q1Δx/2），f（x0+q1Δx/2），…，f（x0+qNΔx/2）可给出函数f（x）在点x0处的一阶导数估值。由泰勒展开有：

图1 估算一阶导数的对称计算节点示意图

油气成藏理论与勘探开发技术

油气成藏理论与勘探开发技术

两式相减，得

油气成藏理论与勘探开发技术

略去（1）式中误差项O（（qiΔx/2）2N+1），将精确的f（i）（x0）换成f（i）（x0）的估计 （x0），并写成矩阵形式，则有

油气成藏理论与勘探开发技术

记

油气成藏理论与勘探开发技术

油气成藏理论与勘探开发技术

设A-1为A 的逆，即 AA-1=I，则（A-1）TAT=I，即（A-1）T为 AT的逆，从而有 AT（A-1）T=I，右乘向量e1得

AT（A-1）Te1=e1 （3）

（3）由方程（2）知

油气成藏理论与勘探开发技术

记向量（A-1）Te1=（c1，c2，…，cN）T，代入式（4）则有

油气成藏理论与勘探开发技术

且由式（3）及（A-1）Te1=（c1，c2，…，cN）T知cn（n=1，2，…，N）满足方程

油气成藏理论与勘探开发技术

式（5）即为用对称任意间距节点计算一阶导数的计算式，cn（n=1，2，…N）通过求解方程（6）得到。若取qn=2n-1（n=1，2，…，N），则cn（n=1，2，…，N）即为等距节点情况的一阶导数具有2N阶精度差分近似的差分系数。

将f（x）做傅氏变换变换到波数域，即f（x） F（k），对f（x）一阶导数进行傅氏变换则有f（1）（x） .i2πkF（k），即f（x）可看作是对f（x）进行的空间线性滤（1）波结果，其滤波响应为i2πk。用f（x）离散点的值估算∂f/∂x得 （x），亦看作是对f（x）进行的空间线性滤波结果，记其滤波响应为 （k），由式（5）知 （k）= cnsin（πkqnΔx）。记

油气成藏理论与勘探开发技术

显然，用 （x）作为f（1）（x）的估计，其精度可由 （k）逼近i2πk的程度，即e（k）来反映。

2 一阶导数的变网格有限差分估算

采用变网格计算，以x方向的一阶导数估算为例来说明两种网格过渡处的处理方法。设细网格节点间距为Δx，粗网格节点间距为mΔx，即粗网格节点间距是细网格节点间距的m倍。为表述方便，不妨取m=3，采用10阶精度，即N=5，如图2所示。图中 表示粗网格节点； 表示粗网格中的加密节点，其与粗网格节点一起组成细网格节点；“+”为求导点位置。①～⑥为粗细网格过渡区的6个一阶导数计算点，计算该6个点处的一阶导数所用的网格节点分别在图2中标示，其中①及其左边大网格区的各计算点按粗网格计算；⑥及其右边的细网格区各计算点按细网格计算；②～⑤4个点需按式（5）计算。

图2 两种网格过渡区计算示意图

一般取2N阶差分计算精度时，在粗、细网格分界的细网格一侧共有N-1个计算点需要采用如式（5）给出的对称不等间距节点差分算子进行一阶导数估算，每一个点的差分算子不同，加上粗网格和细网格内部计算点的差分算子，这样共有N+1个不同的差分算子。N+1个不同的差分算子可事先算好存于内存数组中供每次计算调用。

图3 10阶精度（N=5）不同计算点差分算子对应的滤波响应

①～⑥分别为图2中对应计算点的差分算子的滤波响应；⑦为理想的滤波响应

图3 给出的是图2中的①～⑥求导计算点的估算一阶导数差分算子的滤波响应。对于计算点⑥，由于参加计算的节点为细网格中的相邻2N个网格节点，显然其差分算子的滤波响应与理想的滤波响应最接近，精度最高；对于计算点①，参加计算的节点为粗网格中的相邻2N个网格节点，其差分算子的滤波响应仅在低波数与理想的滤波响应接近，精度最低；而对于计算点②～⑤，精度介于计算点①与⑥之间，且随着参加计算的近距离点增加，其精度越来越高。

在粗细网格过渡区一阶导数的估计精度不一致可能会引起反射，附加计算噪声，现对此做出分析。

用有限差分方法估算空间一阶导数的精度由差分频散决定，有限差分算子的滤波响应只能在低波数段逼近一阶导数算子的滤波响应，在高波数段将出现严重的频散。设kα为具有可接受差分频散水平的最高波数，即在k≤kα的低波数区， （k）与i2πk足够逼近，频散很小，对最终正演结果的影响可忽略不计。记kα=αkN，这里kN=1/（2Δx），0＜α＜1，α值由差分算子决定，如：对均匀网格时的10 阶精度的差分算子，可认为α约为0.6。又设模拟地震波场的最大波数为kmax，它由子波最高频率fmax与最低速度vmin决定，即kmax=fmax/vmin。为了使有限差分数值模拟计算结果不受差分频散的显著影响，必须选取Δx满足α/（2Δx）=kα＞kmax=fmax/vmin，即

油气成藏理论与勘探开发技术

式（7）是保证数值模拟计算精度的空间离散步长Δx应满足的条件，只要满足条件（7），则在波数范围［0，kmax］内差分频散误差可忽略。

采用变网格计算，粗网格尺度也应该满足式（7），即

油气成藏理论与勘探开发技术

否则波场在粗网格区传播就会产生强的频散噪声。式（8）成立，则条件（7）必然满足，由前面对两种网格过渡区的差分算子精度的分析结论知，此时过渡区的差分算子均满足精度要求，即kα＞kmax。既然所有差分算子都能够满足不产生明显的差分频散噪声的要求，差分算子变化也将不会引起明显的计算误差。

因此，对于两种网格过渡区波场传播的衔接问题，上述分析的结论是只要粗网格计算满足精度要求，即不产生明显的差分频散，则按本文的计算方法，在两种网格的过渡区不会出现明显的反射噪声，波场传播的模拟精度基本不受计算网格变化的影响。

3 d性波方程变网格有限差分正演模拟的实现

3.1 实现方法

油气成藏理论与勘探开发技术

以二维d性波方程数值模拟计算来说明变网格有限差分法的实现。利用二维d性波方程进行有限差分正演计算，在交错网格中的递推算式为

油气成藏理论与勘探开发技术

上式中， 和 分别表示 （x，z，t） 和 ∂z（x，z，t） 的有限差分估计，余类推。在式（10）中，ρ，c11，c13，c33及c55都是空变的，为书写简洁，略去了其空间下标。计算中对边界的处理是，在顶边界采用自由边界条件，在其他数值边界则采用PML吸收边界条件［4，6］。

图4 粗、细两种网格中计算节点的布排

在一定的网格尺度下，交错网格计算具有精度高的优点。当采用变网格，即在模型局部用细网格剖分进行计算时，关键是要合理布置计算节点，以保证在所有计算节点均能采用交错网格计算方式。图4 以粗细网格尺度比等于3 为例，给出本文采用的计算节点布排方式，不难看出，要保持交错网格的计算方式，x和z方向粗细网格尺度比 mx和 mz都必须为奇数。

对波场递推计算中空间一阶导数的有限差分估算，根据节点所处位置，分为 3种情况分别处理（为简化图示，不妨设差分算子长度N=2，即4阶精度）：

（1）对细网格区域之外（图4实线框之外）的计算节点采用粗网格进行计算；

（2）对细网格区域内部（图4虚线框内）的计算节点采用细网格进行计算；

（3）对细网格边界区域（图4虚线框与实线框之间）的计算节点按图2描述的方法进行计算，但在计算前对参与计算的某些点（图4 中空心点所在的位置）的波场值需要通过内插得到，本文采用拉格朗日插值方法，内插阶数与差分算子的阶数相同。

3.2 算例

首先，用一个均匀介质模型检验网格变化是否对模拟波场传播产生明显的影响。

如图 5 所示，模型大小为 1000m×1000m，vP=3000m/s，vS=1800m/s，ρ=2.4g/cm3，vz震源位于模型正中位置，采用两种网格划分进行正演计算。第一种网格划分如图5（a）所示，在225～275m深度段（图中阴影部分）用1m×1m的细网格，其余部分用5m×5m的粗网格，采用该变网格进行正演计算，图5（b）和（c）分别是160ms时刻的x分量和z分量波场切片；第二种网格划分如图5（d）所示，在725～775m深度段（图中阴影部分）用1m×1m的细网格，其余部分用5m×5m的粗网格，采用该变网格进行正演计算，图5（e）和（f）分别是160ms时刻的x分量和z分量波场切片。比较两组计算结果可以看到，图（b）与（e），（c）与（f）几乎没有差别，即局部采用细网格对波场传播没有带来明显的影响。为了看得更清楚，将沿过震源的垂直线上（图5（d）中虚线所示）的波传播情况显示于图6，从向上、向下两个方向传播的波的一致性可知，波通过细网格区域时没有产生反射噪声。

图5 变网格计算与波场切片

其次，通过一个含不同尺度圆形洞的介质模型正演结果观察变网格差分计算的效果。

如图7所示，模型大小为1500m×600m，vP=5000m/s，vS=3000m/s，ρ=2.6g/cm3，在模型中部500m深度上分布3个半径分别为5m，10m，20m的3个圆形洞，洞间相距200m，洞内vP=1800m/s，ρ=1.2g/cm3。模拟野外地震观测，炮间距10m，道间距5m，排列长度为1500m，中间激发，激发震源处于模型中部地表，激发子波主频为40Hz，采用变网格计算，粗网格尺度为5m×5m，细网格（图7上的白色虚线框内）尺度为1m×1m，共得150个炮记录，抽出零偏移距剖面，其结果如图8所示。图9（a）是激发点位于模型中部的一炮记录；图9（b）是与图9（a）相对应的对整个模型用1m×1m细网格计算的炮记录。将这两个炮记录进行对比可看到二者差别甚微，证实采用变网格计算基本没有降低计算精度。但采用变网格计算却极大地减少了计算量，就该模型来说，如果对整个模型采用1m×1m的网格进行计算，其计算量约为采用上述变网格计算的计算量的15倍。不难看出，如果对更大尺度的模型进行正演计算，采用变网格计算对减少计算量就显得更有必要。

图6 变网格对波场传播影响测试

图7 含圆洞的介质模型

4 结论

在对复杂介质进行数值模拟计算时，采用变网格计算是解决计算精度与计算量之间矛盾的有效途径。理论分析和数值模拟结果表明，本文给出的适应变网格计算的d性波方程正演模拟实现方法保持了交错网格高阶差分的精度高特点，在粗细网格过渡区不产生人为反射噪声，是一种在保证计算精度前提下减少计算量的有效实用方法。

图8 圆洞模型的零偏移距剖面

图9 变网格与固定小网格计算的炮记录对比

参考文献

［1］董良国.一阶d性波方程交错网格高阶差分解.地球物理学报，2000，43（3）：411～419.

［2］牟永光，裴正林.三维复杂介质地震数值模拟.北京：石油工业出版社，2005.

［3］Jastram C，Tessmer E.Elastic modeling on a grid with vertically varying spacing.Geophys.Prosp.，1994，42：357～370.

［4］Francis Collino，Chrysoula Tsogka.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in an isotopic heterogeneous media.Geophysics，2001，66（1）：294～307

［5］Cunha C A.Elastic modeling in discontinuous media.Geophysics，1993，58（12）：1840～1851.

［6］Rune Mittet.Free-surface boundary conditions for elastic staggered-grid modeling schemes.Geophysics，2002，67（5）：1616～1623.

要在SolidWorks中制作围栏图纸并填充方形网格，可以按照以下步骤 *** 作：

1. 打开SolidWorks并打开你想要制作围栏图纸的图纸文件。

2. 将绘图区域切换到草图模式。如果当前还没有草图，可以通过选择草图工具启动一个新的草图。

3. 选择矩形工具（或者矩形角度工具），绘制一个符合围栏形状的矩形。

4. 在矩形的属性栏中，将线条颜色和线型设置为所需的值，然后手动在矩形内绘制一个方形网格。

5. 选择矩形填充工具，填充整个矩形区域。你可以在矩形填充工具属性栏中选择所需的颜色，并使用“交错”选项创建方形网格。

6. 将绘图区域切换回草图模式，将矩形和内部方形网格复制到正式的围栏图纸中。

7. 如果需要修改网格的大小或形状，可以编辑草图并再次使用矩形工具进行绘制。填充可用同样演示7中的方法。

希望这些步骤能够确保你能够成功地创建围栏图纸并填充方形网格。

---