复杂构造网格化及高精度地震波场谱元法数值模拟

李洪建，韩立国，巩向博

(吉林大学地球探测科学与技术学院，吉林长春130026)

随着油气勘探开发的不断深入，油气藏勘探逐渐从常规构造油气藏向复杂隐伏构造油气藏和隐蔽油气藏发展，复杂的地下地质结构和地层条件对地震勘探方法技术提出了更高的要求。地震波场数值模拟是研究复杂构造地震波传播问题的重要工具，人们充分利用现有的地质和地球物理资料，不断提高数值模拟的精度，为油气藏勘探开发提供更准确的依据[1]。

正演数值模拟在油气勘探中应用广泛，李军锋[2]分析了高精度有限差分波场正演模拟方法，并对张性断层和逆掩推覆断层进行了模拟计算；王秀明等[3]利用高阶交错网格有限差分法模拟了地震波在烃类储集层中的传播；刘春园等[4]利用变网格有限差分法对孔洞模型等碳酸盐岩储层进行了研究。有限差分法计算速度快，但网格剖分方法无法适应不规则异常地质界面的情况，且存在数值频散。刘炯等[5]提出了错格伪谱法并模拟了不均匀缝洞型储层介质中弹性波的传播问题；刘渊博等[6]提出了分段点最佳化的网格优化算法，提高了伪谱法对非光滑最优控制问题的求解效率。伪谱法模拟精度高，但伪谱法网格剖分同样无法很好地适应边界的变化。刘广峰等[7]利用有限元方法开展了油气储层地应力模拟研究；王月英[8]采用MPI并行计算方法，提升了有限元正演模拟数据容量和计算速度。有限元法可以很好地模拟复杂边界构造，但大规模的计算量使它难以发挥优势。

谱元法数值模拟方法结合了伪谱法高精度和有限元法灵活处理复杂地质结构的优点，在模拟复杂构造方面具有很大的优势，在应用到地震勘探领域的20多年里，国内外专家、学者对谱元法数值模拟技术及其应用做了一系列的改进，使其更具有实用性。Patera[9]最早将谱元法应用于流体动力学计算；Komatitsch等[10-17]编写了从一维到三维的谱元法计算程序，在地震波场计算中取得了很好的应用效果；王秀明等[18]提出了逐元算法，大大提高了谱元法计算效率。但是，常规的谱元网格剖分方法仍受到地层厚度限制，导致部分区域网格剖分过于致密或稀疏，模拟精度受到很大影响。对于复杂构造的地震正演模拟，发展一种既快速方便，又能灵活适应地层边界变化的正演算法尤为重要。

在谱元法数值模拟中，灵活使用四边形网格对模型进行空间离散，能更好地逼近复杂构造。在三角形网格内部插入一个节点，连接这一节点和三角形质心，网格可以被剖分为3个四边形，但这种方法剖分出的四边形往往形状不规则，单元质量不高；对于合适的几何域，四边形映射网格剖分(quadrilateral mapped meshing)[19]或子映射(sub-mapping)[20]方法可以将整个区域离散为一组结构化的四边形，每个内部节点都连接有4个单元，这种方法快速且有效，但是仅仅适用于某些特定的问题。为此，更为普遍的非结构化算法被提出。非结构化四边形网格算法可以分为直接法和间接法两类。直接法四边形网格剖分方法可以借助于域分解[21-22]把区域直接剖分为简单的可变换凸四边形单元，但是其剖分结果存在大量不稳定单元。区别于直接法，间接法四边形网格剖分方法首先将域剖分为三角形网格，进而合成四边形网格。Johnston等[23]使用拆分和交换方法，三角形单元被尽可能全部重组为四边形单元，但仍有一部分三角形网格余留；Lee等[24]和Lo[25]提出了首先在边界处将三角形重组形成四边形单元，进而利用波前法在区域内部依次将剩余所有的三角形重组为四边形单元的方法，这种重组方法得到的四边形网格质量较高，剖分效果较好。

我们研究并提出了复杂构造网格剖分及其高精度谱元数值模拟方法。首先基于Delaunay三角网格剖分方法将区域剖分为三角形网格；然后采用波前法将所有三角形网格重组，得到谱元法所需的四边形网格。应用基于复杂构造网格剖分的高精度谱元法对Marmousi速度模型和二维随机介质模型进行数值模拟，并与常规谱元法的模拟结果进行对比，以验证其有效性和应用效果。

1 方法原理

根据Galerkin加权余量原理，与微分方程初、边值问题对应的波动方程弱形式可以表示为

(1)

将模型剖分为若干个离散的单元，在每个单元上做离散化运算。对于测试函数空间内的位移值uh，有

(2)

式中：gh表示边界位移；vh表示除边界外测试函数空间的位移值。用截断的Chebyshev正交多项式或者Legendre多项式确定单元基函数NA(x)，测试函数wh，vh，gh在单元内各个节点上的节点值分别为ciA(t)，diB(t)，giB(t)，则有

(3)

式中：η为全局节点集合；ηgi为满足边界条件的节点集合。经过计算，谱元方程矩阵形式可以简写为

(4)

(5)

式中：ei，ej为单位向量；fi为单元体积的体力;hi为边界应力；cijkl为弹性系数。

图1 Delaunay三角形网格重组示意图解

针对某一起伏地表层状模型，图2给出了常规谱元法网格剖分结果和本文复杂构造网格剖分法的剖分结果。从图2中可以看出，在起伏地表和起伏界面处，二者都能很好地适应边界的变化；但对比可见图2a剖分效果相对较差，因为常规谱元法的网格剖分方法必须保证每层具有特定的网格数目，由于地层厚度的变化，难免会导致部分区域网格剖分过于致密或稀疏，从而导致网格大小不一，差别很大；而图2b则有效地避免了这一问题的产生，四边形网格大小均匀，质量较好。

图2 起伏地表层状模型的常规谱元法(a)和复杂构造网格剖分方法(b)的网格剖分结果

Delaunay三角网格剖分的任意性导致重组之后的网格形状为不规则四边形，无法直接用于运算，需要通过参数变换，才能将单元运算法则应用到形状不规则的网格单元中。将参考单元映射到实际物理区域的变换方法称为等参变换，即单元几何形状变换，单元内场函数采用相同的节点参数及相同的插值函数，参与变换的单元称为等参单元。图3为具有4个节点单元的等参变换的基本示意图。

图3 4个节点单元的等参变换示意图解

为将自然坐标(ξ，η)中形状规则的单元转换成笛卡尔坐标(x,y)中形状扭曲的单元，需要建立一个坐标变换：

(6)

并建立以下映射：

(7)

式中：Na表示基函数；nen表示单元中的节点数，此例中nen=4。由于

(8)

式中：Na,ξ，Na,η分别表示基函数Na在ξ，η方向的导数。通过推导，得到基函数Na关于x，y的导数：

(9)

Na,ξ和Na,η来自于几何形状规则的单元计算,容易求得。并且

(10)

其中，J为Jacobian矩阵。

(11)

即可得

(12)

通过等参变换，即可以利用规则坐标来计算我们剖分所得的不规则网格。把所有单元按照一定的顺序进行组装，求解总体谱元方程，即能得到全局的近似解。

2 数值模拟算例

2.1 Marmousi速度模型

首先利用基于复杂构造网格剖分方法的高精度谱元法和常规谱元法对Marmousi速度模型(图4a) 进行模拟对比研究，模型大小为4170m×1810m。为了测试复杂构造网格剖分方法的优势，按照大致的速度边界将整个模型区域划分为15个区域，在每个速度区域进行Delaunay三角网格剖分和间接法重组，得到谱元法计算所需要的四边形网格(图4b)，网格总数约为120000个。由于网格剖分较细，将图4b中黑框区域放大如图5a，再将图5a中黑框区域放大如图5b。对比图6 所示的常规谱元法剖分结果(将模型区域剖分为400×300=120000个网格)，可以清楚地看出，复杂构造网格剖分方法不但能很好地适应复杂构造速度边界和倾角的变化，且网格均匀，质量较高，对复杂地层适应性更好。

图4 Marmousi速度模型(a)及其复杂构造网格剖分法的剖分结果(b)

图5 Marmousi模型复杂构造网格剖分法剖分结果(图4a)按黑框区域的依次放大显示a 图4b中黑框区域的放大显示; b 图5a中黑框区域的放大显示

图6 Marmousi模型的常规谱元法网格剖分结果(a)及其黑框区域的放大显示(b)

选取频率为40Hz的雷克子波作为震源，激发点位置在模型表面(2085m,25m)处；在(100m,25m)到(4100m,25m)处均匀排列共251个接收道，道间距为16m；取采样点数9000，采样率0.2ms，总记录长度为1.8s，分别采用基于复杂构造网格剖分方法的高精度谱元法和常规谱元法对Marmousi速度模型做正演数值模拟。图7为高精度谱元法数值模拟0.5，0.6，0.7，0.8s时刻的波场快照；图8 为高精度谱元法和常规谱元法正演模拟的单炮记录对比结果；图9是对应图8中黑框区域(50～240道，4500～9000采样点(0.9～1.8s))的放大显示对比。

从高精度谱元法数值模拟的波场快照(图7)可以看出，弹性波遇到Marmousi模型的速度边界时，在界面处发生较清晰的反射和透射，波形清晰连续。但由于该模型比较复杂，其反射、透射也很复杂。

从模拟单炮记录(图8)及其放大结果(图9)可明显看出，与常规谱元法相比，基于复杂构造网格剖分的高精度谱元法由于网格剖分更适应实际速度界面，同相轴更清晰、连续，且频散很小，模拟精度高，具有很高的信噪比和分辨率。

图7 Marmousi模型高精度谱元法数值模拟不同时刻的波场快照a 0.5s; b 0.6s; c 0.7s; d 0.8s

图8 Marmousi模型的高精度谱元法(a)和常规谱元法(b)正演模拟单炮记录

图9 Marmousi模型的高精度谱元法(a)和常规谱元法(b)正演模拟单炮记录中黑框区域的放大显示

2.2 随机介质模型

在石油勘探中，实际储层(如碳酸盐岩储层)往往是非均匀的，大量微小异常极不规则地分布在非均匀介质中，正演模拟中如果采用波长远大于异常尺度的地震波，非均匀异常对地震波的传播影响较小。但在如今高精度的地球物理探测中，对于高、宽频地震信号，这些小尺度异常却无法忽略，因为它们对地震波的传播产生十分明显的影响。为了验证复杂构造高精度谱元法的应用效果，建立如图10 所示的二维水平地表层状随机介质模型。模型大小为1200m×1200m，其中随机介质的物理参数如表1所示。

图10 层状随机介质速度模型

采用基于复杂构造网格剖分的高精度谱元法将层状随机介质模型区域剖分为大约90000个网格，结果示于图11a，由于网格剖分较细，图11b给出了图11a中心区域的局部放大显示。对比图12所示的常规谱元法网格剖分结果(将模型区域剖分为300×300=90000个网格)，可以看出，复杂构造网格剖分法的剖分结果大小均匀，质量较高，对复杂界面适应性更好。

为考虑非均匀性对地震波的影响，地震波长要尽可能接近异常的尺度，随机介质模型中非均匀异常不规则且尺度相对微小，因此需要选用高频地震信号进行谱元法数值模拟。在模型中心(600m,480m)位置处选取频率为150Hz的雷克子波作为震源，并在(250m,480m)到(950m,480m)均匀排列201个检波器，取采样点数3360，采样率0.1ms，得到每个时刻的波场快照和单炮记录。图13为高精度谱元法数值模拟0.08，0.12，0.16，0.20s时刻的波场快照；图14为二维随机介质速度模型高精度谱元法和常规谱元法正演模拟单炮记录对比；图15为对应图14中黑框区域(25～201道，2400～3360采样点(0.240～0.336s))的放大显示。

从高精度谱元法数值模拟的波场快照(图13)

中可以看出，直达波的波前面不再是在均匀介质中传播时呈现的圆形，可见随机介质中存在速度方面的微小差异，具有明显的不均匀性；在速度模型内部明显存在各种微小的散射波，它们由非均匀地质体的速度差异和密度差异引起。

表1 随机介质物性参数

图11 随机介质模型的复杂构造网格剖分法剖分结果(a)及其中心区域的局部放大显示(b)

图12 随机介质模型的常规谱元法网格剖分结果(a)及其中心区域的局部放大显示(b)

图13 随机介质模型高精度谱元法数值模拟不同时刻的波场快照a 0.08s; b 0.12s; c 0.16s; d 0.20s

图14 随机介质模型的高精度谱元法(a)和常规谱元法(b)正演模拟单炮记录

图15 随机介质模型高精度谱元法(a)和常规谱元法(b)模拟单炮记录中黑框区域的放大显示

从图14和图15的模拟单炮记录来看，与常规谱元法相比，基于复杂构造网格剖分的高精度谱元法模拟结果具有更高的信噪比和分辨率，同相轴连续性好，精度更高。高精度谱元法能将微观层次的有效地震信息更清晰地反映在地震记录剖面上，有助于在高精度地震勘探中发现小尺度的储层非均质体。

3 结论和讨论

在空间上基于Delaunay三角网格剖分方法对模型区域进行剖分，再使用间接法首先在边界处将三角形重组，进而利用波前法在区域内部依次将剩余所有的三角形重组为谱元法计算所需的四边形单元，形成了适应复杂构造模拟网格剖分的新方法。与常规网格剖分方法相比，复杂构造网格剖分方法不但能灵活地适应速度边界的变化，且网格大小均匀、质量较高，有效地解决了常规谱元法数值模拟在复杂构造区域存在网格剖分过密或过疏的不足。

对Marmousi模型进行谱元法数值模拟的研究结果表明，在复杂网格空间离散条件下，基于复杂构造网格剖分的高精度谱元法能够更好地模拟复杂构造，有效信号清晰且连续，数值频散小，精度很高。随机介质的物理性质(如速度、密度等)在微观层次的微小差异导致介质内部存在波阻抗差，并影响其散射波场的强弱和分布，针对二维层状随机介质模型的高精度谱元法数值模拟可以将这些有效信息更清晰地记录下来。数值模拟结果展示了基于复杂构造网格剖分的高精度谱元法在高精度地震勘探中具有较好的应用前景，提高了勘探精度，有助于发现小尺度的储层非均质体。

参 考 文 献

[1] 王锡文,秦广胜,赵卫锋,等.正演模拟技术在地震采集设计中的应用[J].地球物理学进展,2012,27(2)：642-650

Wang X W,Qing G S,Zhao W F,et al.Application of forward modeling in seismic acquisition design[J].Progress in Geophysics,2012,27(2)：642-650

[2] 李军峰.高精度有限差分波场正演模拟方法[D].成都：成都理工大学,2012

Li J F.High precision wave field forward modeling of finite difference method[D].Chengdu：Chengdu University of Technology,2012

[3] 王秀明,张海澜,王东,等.利用高阶交错网格有限差分法模拟地震波在非均匀孔隙介质中的传播[J].地球物理学报,2003,46(6)：842-849

Wang X M,Zhang H L,Wang D,et al.Simulations of seismic wave propagation in inhomogeneous porous media using high-order staggered finite difference method[J].Chinese Journal of Geophysics,2003,46(6)：842-849

[4] 刘春园,朱生旺,魏修成,等.随机介质地震波正演模拟在碳酸盐岩储层预测中的应用[J].石油物探,2010,49(2):133-139

Liu C Y,Zhu S W,Wei X C,et al.Random medium seismic forward modeling in carbonate reservoir prediction[J].Geophysical Prospecting for Petroleum,2010,49(2):133-139

[5] 刘炯,马坚伟,杨慧珠,等.缝洞型储层错格伪谱法模拟[J].石油地球物理勘探,2008,43(6):723-727

Liu J,Ma J W,Yang H Z,et al.Vuggy reservoir simulation staggered grid pseudo-spectral method[J].Oil Geophysical Prospecting,2008,43(6):723-727

[6] 刘渊博,朱恒伟.伪谱法求解非光滑最优控制问题的网格优化[J].系统工程与电子技术,2013,35(11)：2396-2399

Liu Y B,Zhu H W.Pseudo-spectral method for solving non-smooth mesh optimization optimal control problems[J].Systems Engineering and Electronics,2013,35(11)：2396-2399

[7] 刘广峰,陆红军,何顺利.有限元法开展油气储层地应力研究综述[J].科学技术与工程,2009,9(24):7430-7435

Liu G F,Lu H J,He S L.A Summary of oil and gas reservoirs to stress finite element method[J].Science Technology and Engineering,2009,9(24):7430-7435

[8] 王月英.基于MPI的三维波动方程有限元法并行正演模拟[J].石油物探,2009,48(3)：221-243

Wang Y Y.3-D Finite dimensional wave equation forward simulation based on MPI parallelmodeling[J].Geophysical Prospecting for Petroleum,2009,48(3)：221-243

[9] Patera A T.A spectral element method for fluid dynamics:laminar flow in a channel expansion[J].Journal of Computational Acoustics,1994,2(4)：371-422

[10] Komatitsch D.Spectral and spectral-element methods for the 2D and 3D elasto-dynamics equations in heterogeneous media[D].Paris:Insti-tut de Physique du Globe,1997

[11] Komatitsch D,Vilotte J P.The spectral element method:an efficient tool to simulate the seismic response of 2D and 3D geological structures[J].Bullein of the Seismolgical Society of America,1998,88(2):368-392

[12] Tromp J,Komatitsch D.The spectral element method:spectral-element and adjoint methods in seismology[J].Communications in Computational Physics,2008,3(1):1-32

[13] Lee S J,Komatitsch D.Effects of topography on seismic-wave propagation:an example from Northern Taiwan[J].Bulletin of the Seismological Society of America,2009,99(1):314-325

[14] Lee S J,Chan Y C,Komatitsch D.Effects of realistic surface topography on seismic ground motion in the Yangminshan region of Taiwan based upon the spectral-element method and LiDAR DTM[J].Bulletin of the Seismological Society of America,2009,99(2):681-693

[15] Komatitsch D,Erlebacher G,Göddeke G,et al.High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster[J].Journal of Computational Physics,2010,229(20):7692-7714

[16] Komatitsch D,Tromp J.Introduction to the spectral element method for three-dimensional seismic wave propagation[J].Geophysical Journal International,1999,139(3):806-822

[17] Komatitsch D,Tromp J.Spectral-element simulations of wave propagation in porous media[J].Geophysical Journal International,2008,175(1):301-345

[18] 王秀明,Serianig G,林伟军.利用谱元法计算弹性波场的若干理论问题[J].中国科学(G辑:物理学 力学 天文学),2007,37(1)：41-59

Wang X M,Serianig G,Lin W J.Using spectral element method to calculate some theory problems about elastic wave field[J].Science in China(Physics,Mechanics & Astronomy),2007,37(1)：41-59

[19] Cook W A,Oakes W R.Mapping methods for generating three-dimensional meshes[J].Computers in Mechanical Engineering,1982,67-72

[20] White D R.Automated hexahedral mesh generation by virtual decomposition[J].Meshing Roundtable,1985,165-176

[21] Baehmann P L,Wittchen S L,Shephard M S,et al.

Robust geometrically-based,automatic two-dimensional mesh generation[J].International Journal for Numerical Methods in Engineering,1987,24:1043-1078

[22] Talbert J A,Parkinson A R.Development of an automatic,two dimensional nite element mesh generator using quadrilateral elements and Bezier curve boundary dentition[J].International Journal for Numerical Methods in Engineering,1991,29:1551-1567

[23] Johnston B P,Sullivan J M,Kwasnik A.Automatic

conversion of triangularnite element meshes to quadrilateral elements[J].International Journal for Numerical Methods in Engineering,1991,31:67-84

[24] Lee C K,Lo S H.A new scheme for the generation of a graded quadrilateral mesh[J].Computers and Structures,1994,52(5):847-857

[25] Lo S H.A new mesh generation scheme for arbitrary planar domains[J].International Journal for Numerical Methods in Engineering,1985,21:1403-1426