基于双二次插值的地震波场有限元法数值模拟

周建科，印兴耀，曹丹平

(中国石油大学(华东)地球科学与技术学院，山东青岛266580)

波动方程正演模拟对认识地震波传播规律和指导地震资料解释等具有重要意义[1-2]，要想精确地模拟地震波在地下介质中的传播，一方面要采用较为接近实际地层的地球物理模型，另一方面要采用高精度的数值模拟算法。有限差分法(FDM)和有限元法(FEM)是地震波场正演模拟中非常重要的两种方法，有限差分法具有编程简单，计算效率高等特点，在地震模拟中得到广泛的研究和应用，但难以处理起伏地表和自由边界条件[3-4];有限元法能够适应剧烈的物性变化，自然地满足自由边界条件，因此具有精度高、可模拟任意复杂结构、易于进行边界处理等优点[5-8]，其缺点是计算量大和内存占用高。

有限元法具有多种网格剖分方式，能够对任意复杂的边界进行有效剖分，因而能够准确模拟起伏地形、复杂构造和复杂介质条件下的地震波场，种类繁多的插值函数可以提供不同精度的数值计算结果。当单元采用线性函数进行插值时，得到的数值结果往往精度较低，达不到预期目标，当采用高次函数进行插值时，具有较高的数值模拟精度，但计算量以及内存的占用也会增加。史明娟等[9]以及刘云等[10]采用双二次插值有限元法(以下简称双二次插值法)实现了大地电磁的数值模拟，得到高精度的数值结果;戴前伟等[11]实现了双二次插值的探地雷达有限元数值模拟，与双线性插值有限元法(以下简称双线性插值法)相比，双二次插值法的数值模拟结果更能突出异常体的响应。

应用有限元法解波传播问题时需占用巨量的机时和内存，使得它在地球物理中的应用受到极大限制。郭建[12]和刘有山等[8]根据结构刚度矩阵稀疏性的特点，采用稀疏矩阵的压缩存储行(Compressed Sparse Row，CSR)格式。在此不仅考虑了结构刚度矩阵的稀疏性，还考虑其对称性，采用紧凑存储法只存储结构刚度矩阵下三角部分的非零元素，同时，在波场的递推过程中，只有结构刚度矩阵每一行的非零元素与对应的节点位移进行相乘;根据质量守恒原则，采用对角的集中质量矩阵代替一致质量矩阵，避免显式有限元法在每个时间步长上进行大型稀疏矩阵的求逆运算。

我们在矩形网格剖分情况下，采用双二次插值法实现了二维声波方程的数值模拟，采用对角的集中质量矩阵和紧凑存储法提高计算效率、减少内存的占用，为实现有限元法的高精度数值模拟提供一定的指导。

1 方法原理

1.1 有限元法求解二维声波方程

使用伽勒金法求解二维标量波动方程得到有限元方程组[13]

(1)

式中：U(t)代表t时刻结构节点位移列向量;M和K分别代表结构的质量矩阵和刚度矩阵：

(2)

式中：Me和Ke分别为单元的质量矩阵和刚度矩阵;N为结构中单元的个数;c代表单元内介质的声波速度;N为形函数矩阵;∑代表单元质量(刚度)矩阵组装为结构质量(刚度)矩阵的法则。

对于(1)式，时间因子采用二阶中心差分格式进行求解，并考虑震源的加载，得到

(3)

式中：F(t)为t时刻作用在节点上的等效载荷;Δt为时间采样间隔。

进行波场的递推时，需要给定初始条件：

(4)

1.2 单元内双二次插值

首先将求解区域Ω剖分成矩形单元，在单元内进行双二次插值。在每个单元上取8个节点：4个角点和4条边的中点，8个节点的编号次序及坐标如图1所示。母单元与子单元间的坐标变换关系为

(5)

其中，x0,y0是子单元中点的坐标，a与b分别为子单元的长和宽。

采用双二次插值，母单元内任意一点的位移u可以表示为

(6)

式中：a1,…,a8为待定常数;ξ,η为母单元上的坐标。

将母单元8个节点的坐标以及相应的位移代入(6)式，整理可得

(7)

图1 单元剖分示意a 子单元; b 母单元

(8)

将(5)式、(8)式代入(2)式，采用等参单元和高斯数值积分[14]，得到单元刚度矩阵Ke和一致质量矩阵Me。

两矩阵计算式分别为

1.3 集中质量矩阵

在地震波场数值模拟中，一般采用对角的集中质量矩阵代替一致质量矩阵，以避免时间循环过程中的求逆运算。对线性插值单元集中质量矩阵的求取，一般将一致质量矩阵每一行所有元素叠加到对角线上，然后令非对角线元素为零。但对高次插值单元而言，这样的做法不可取，例如将(10)式中的矩阵的每一行元素叠加到对角线上后，4个角节点分配到的质量为负值，这在力学上是不合理的，同时也将对计算精度产生不良影响[15]。

对于矩形双二次插值单元，采用质量守恒原理来确定单元集中质量矩阵。设每个角节点分配到的单元质量为m1，每条边的中点分配到的单元质量为m2，单元的总质量为m。定义质量分布参数β满足

(11)

根据质量守恒定理，有

(12)

将(12)式代入(11)式整理得到：

(13)

(14)

矩形单元双二次插值的集中质量矩阵的形成不是唯一的，例如以单元一致质量矩阵的对角线元素为权重时有β=3/76=0.0395，此外取β=1/36=0.0278在实际中也运用得较多[15]。本文进行数值模拟时，β的取值为3/76，即β=0.0395。

1.4 结构刚度矩阵的紧凑存储法

根据结构刚度矩阵的集成法则(附录A)可知，结构刚度矩阵每一行只有少量的非零元素，且具有对称性，每个节点的位移只受到与其共点单元节点的影响。图2是计算区域中的任一共点单元，可见，对于非边界角节点(如节点L3)而言，对其位移产生影响的节点只有21个(包括该节点本身)，而对非边界中点的位移产生影响的节点只有13个，与边界节点位移相关的节点则更少。因此由双二次插值法得到的结构刚度矩阵每一行非零元素最多为21个，再根据其对称性，只需存储结构刚度矩阵下三角部分的非零元素即可。根据以上分析，本文采用紧凑存储法来存储结构刚度矩阵，需要3个一维数组：浮点型数组CS按行存储结构刚度矩阵下三角部分的非零元素，整型数组GA存放CS中对应元素的列号，整型数组ID存放下三角部分每一行非零元素的个数。

以水平方向有nx个单元，垂直方向有nz个单元的均匀剖分结构为例，结构节点与单元编号规则为：先从左到右，然后由上至下，单元节点的局部编号如图1所示。下面以节点L3为例来说明紧凑存储法的原理。节点L3决定结构刚度矩阵的第L3行，该行非零元素的列号为与其共点单元节点的结构编号，再根据结构刚度矩阵的对称性，第L3行所需存储的元素仅为11个，其列号分别为：L1-2，L1-1，L1，L1+1，L1+2，L2-1，L2，L2+1，L3-2，L3-1，L3，这11个数用数组GA来储存，以便后续使用。一个具有nx×nz个单元的结构，结构刚度矩阵中元素个数为(3nxnz+2nx+2nz+1)2，而采用紧凑存储法后，结构刚度矩阵中需要存储元素的个数仅为：25nxnz+5nx+5nz+1，辅助数组GA与ID中元素的个数分别为25nxnz+5nx+5nz+1与3nxnz+2nx+2nz+1。

图2 共点单元示意

由(3)式可知，在时间循环过程中，涉及到结构刚度矩阵与结构节点位移列向量的矩阵乘法，如果结构刚度矩阵中的全部元素参与运算，这将非常耗时。众所周知，矩阵中的零元素对矩阵乘法是没有贡献的，利用数组ID和GA可以找到非零元素在结构刚度矩阵中的具体位置，因此在时间循环过程中，只需结构刚度矩阵每一行的非零元素与对应的节点位移进行相乘。变带宽存储法[14]虽然也只存储了结构刚度矩阵下三角部分的元素，但其把半带宽中的零元素也储存了，同时在计算中，半带宽中的零元素也参与了运算，对大型结构不再适用(节点数越多，同一单元节点编号的最大值与最小值之差也就越大，半带越宽，半带宽中的零元素也就越多)。

在内存为4GB，处理器为Intel(R)Core(TM)2 DUO E8400 @3.00GHz的计算机上对比分析了紧凑存储法与变带宽储存法的计算效率以及所占内存，具体结果见表1。当模型网格数为100×100时，变带宽存储法进行300次时间循环的耗时和结构刚度矩阵内存占用分别是紧凑存储法的19.20倍和14.47倍，当网格数为400×400时，相应的倍数变为81.90和56.90。

表1 不同存储方法的CPU耗时以及内存占用统计

2 模型试算

采用C语言编制双二次插值法二维声波方程地震波场数值模拟程序，运行在内存4GB，处理器Intel(R)Core(TM)2 DUO E8400 @3.00GHz的计算机上。震源选用定点上延迟的雷克子波，峰值频率fp=40Hz。

2.1 计算效率比较

为了说明双二次插值法的优缺点(计算效率及内存占用情况)，分别采用双线性插值法和双二次插值法进行地震波场的数值模拟，在无可见数值频散的情况下，比较二者的计算耗时以及存储需求量。在本算例中选取纵波波速c=1900m/s的均匀介质模型，计算区域Ω={(x,z)|0≤x≤2000m,0≤z≤2000m}，时间采样间隔Δt=0.4ms，震源位于模型中央。采用2种算法得到的T=0.4s时刻波场快照如图3所示。其中，图3a至图3c是在不同尺寸单元网格下采用双线性插值法得到的波场快照，当单元网格大小为2.5m×2.5m(图3a)和2.0m×2.0m(图3b)时，波场快照可见微弱的数值频散，当单元网格减小到1.7m×1.7m时(图3c)无可见数值频散，此时结构中的网格数为1176×1176;图3d至图3f则为在不同尺寸单元网格下采用双二次插值法得到的波场快照，在6.0m×6.0m(图3d)以及5.5m×5.5m(图3e)的单元网格下，波场快照具有微弱的数值频散，当单元网格减小到5.0m×5.0m时，数值频散得以消除，此时网格数为400×400。也就是说，当震源主频为40Hz、介质速度为1900m/s时，在保证无可见数值频散情况下，双线性插值法的单元网格最大尺寸为1.7m×1.7m，双二次插值法的单元网格最大尺寸为5.0m×5.0m。

表2给出了本算例在无可见数值频散下，采用紧凑存储法时2种数值模拟方法的CPU耗时(1001次时间循环)以及内存占用情况。可见，为使两者的数值模拟结果均无可见的数值频散，双线性插值法占用的内存是双二次插值法的一倍，且双二次插值法的单步耗时也比双线性插值法少。在这里值得指出的是，双二次插值法的稳定性条件比双线性插值法更为苛刻，因此双线性插值法的时间采样间隔可以取得比双二次插值法大，当两者计算时间相同时，双线性插值法可以选取更大的时间步长，时间循环次数也就变少，有可能使得计算耗时比双二次插值法少，因此本文得出的结论是双二次插值法的单步耗时比双线性插值法少。

表2 无可见数值频散情况下两种数值模拟方法的CPU耗时以及内存占用量

图3 均匀介质模型T=0.4s时刻波场快照a 双线性插值法,单元网格大小2.5m×2.5m; b 双线性插值法,单元网格大小2m×2m; c 双线性插值法,单元网格大小1.7m×1.7m; d 双二次插值法,单元网格大小6.0m×6.0m; e 双二次插值法,单元网格大小5.5m×5.5m; f 双二次插值法,单元网格大小5.0m×5.0m

2.2 复杂构造模型

为了进一步验证双二次插值法对复杂构造的模拟精度，设计了如图4所示的模型，模型包括一个垂直断层和一个背斜构造，长为3600m，深2000m。模型中的阶梯状界面上有4个拐点，其中点1与顶界面以及点2的垂直距离均为400m，点2与点3和4的水平距离分别为400,2000m，点2到其下水平界面的垂直距离为400m，点4与右边界的水平距离为800m，背斜弧高500m，由上至下速度依次为2000，3000，4000m/s。震源位于(1800m，200m)处，四周边界均采用改进的Sarma吸收边界条件[16]，厚度为400m。单元网格大小为4m×4m，考虑到稳定性条件，时间采样间隔取0.2ms，计算时长T=1.6s。

图5给出了0.5，0.6，0.7，0.8，0.9，1.0s时刻的波场快照(单元中心点的波场值由单元8个节点的波场值插值得到)。从图5中可以清楚地看到地震波在4个拐点处产生的绕射波、界面的反射波以及绕射波形成的多次反射波，说明双二次插值法在模拟复杂构造时具有较好的计算精度。单炮记录如图6所示，由于离散背斜界面采用矩形网格，因此可以观察到由阶梯状界面产生的非物理绕射波。

在该数值算例中，实际参与计算(包含边界条件)的网格数为1100×700，进行8001次时间循环耗时4045s，平均单步耗时0.5056s，采用紧凑存储法存储结构刚度矩阵占用的内存为155.76MB，结构集中质量矩阵占用的内存为8.83MB。

图4 复杂构造速度模型

图5 复杂构造模型不同时刻的波场快照a 0.5s; b 0.6s; c 0.7s; d 0.8s; e 0.9s; f 1.0s

图6 复杂构造模型单炮记录

3 结束语

我们采用双二次插值法实现了二维声波方程的有限元法数值模拟，采用质量守恒定律得到对角的集中质量矩阵，结构刚度矩阵采用紧凑存储。模型试算结果表明，双二次插值法能够提高地震波场数值模拟的精度，采用集中质量矩阵和刚度矩阵的紧凑存储后，内存消耗以及计算速度能够满足实际生产要求。双二次插值法和双线性插值法的对比实验表明，在无可见数值频散情况下，双二次插值法消耗的内存更少，单步耗时更短。基于双二次插值法的以上优越性，该算法可望进一步推广到弹性波数值模拟以及偏移成像等领域。

附录A 双二次插值法结构刚度矩阵的组装

根据单元刚度矩阵系数与结构刚度矩阵系数的对应关系(具体对应关系请参考文献[12])，结构刚度矩阵第L3(图2)行需要进行叠加的元素为

参 考 文 献

[1] 吴国忱,王华忠.波场模拟中的数值频散分析与校正策略[J].地球物理学进展,2005,20(1):58-65

Wu G C,Wang H Z.Analysis of numerical dispersion in wave-field simulation[J].Progress in Geophysics,2005,20(1):58-65

[2] 孙林洁,印兴耀.基于PML边界条件的高倍可变网格有限差分数值模拟方法[J].地球物理学报,2011,54(6):1614-1623

Sun L J,Yin X Y.A finite-difference scheme based on PML boundary condition with high power grid step variation[J].Chinese Journal of Geophysics,2011,54(6):1614-1623

[3] Alford R M,Kelly K R,Boore D M.Accuracy of finite-difference modeling of the acoustic wave equation[J].Geophysics,1974,39(6):834-842

[4] 殷文,印兴耀,吴国忱,等.高精度频率域弹性波方程有限差分方法及波场模拟[J].地球物理学报,2006,49(2):561-568

Yin W,Yin X Y,Wu G C,et al.The method of finite difference of high precision elastic wave equations in the frequency domain and wave-field simulation[J].Chinese Journal of Geophysics,2006,49(2):561-568

[5] Marfurt K J.Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations[J]. Geophysics,1984,49(5):533-549

[6] 薛东川,王尚旭,焦淑静.起伏地表复杂介质波动方程有限元数值模拟方法[J].地球物理学进展,2007,22(2):522-529

Xue D C,Wang S X,Jiao S J.Wave equation finite-element modeling including rugged topography and complicated medium[J]. Progress in Geophysics,2007,22(2):522-529

[7] 薛东川，王尚旭.波动方程有限元叠前逆时偏移[J].石油地球物理勘探,2008,43(1):17-21

Xue D C,Wang S X.Wave-equation finite element prestack reverse-time migration[J].Oil Geophysical Prospecting,2008,43(1):17-21

[8] 刘有山,滕吉文,刘少林,等.稀疏存储的显式有限元三角网格地震波数值模拟及其PML吸收边界条件[J].地球物理学报,2013,56(9):3085-3099

Liu Y S，Teng J W，Liu S L，et al.Explicit finite element method with triangle meshes stored by sparse format and its perfectly matched layers absorbing boundary condition[J].Chinese Journal of Geophysics,2013,56(9):3085-3099

[9] 史明娟,徐世浙,刘斌.大地电磁二次函数插值的有限元法正演模拟[J].地球物理学报,1997,40(3):421-430

Shi M J,Xu S Z,Liu B.Finite element method using quadratic element in MT forward modeling[J]. Chinese Journal of Geophysics,1997,40(3):421-430

[10] 刘云,王绪本.大地电磁二维自适应地形有限元正演模拟[J].地震地质,2010,32(3):382-391

Liu Y,Wang X B.FEM using adaptive topography in 2D MT forward modeling[J].Seismology and Geology,2010,32(3):382-391

[11] 戴前伟,王洪华,冯德山,等.基于双二次插值的探地雷达有限元数值模拟[J].地球物理学进展,2012,27(2):736-743

Dai Q W,Wang H H,Fei D S,et al.Finite element numerical simulation for GPR based on quadratic interpolation[J].Progress in Geophysics,2012,27(2):

736-743

[12] 郭建.一种有限元快速算法[J].石油物探,1991,30(2):36-43

Guo J.A kind of fast finite element algorithm[J].Geophysical Prospecting for Petroleum,1991,30(2):36-43

[13] 杜世通.变速不均匀介质中波动方程的有限元法数值解[J].华东石油学院学报,1982,6(2):1-20

Du S T.Finite element numerical solution of wave propagation in non-homogeneous medium with variable velocities[J].Journal of East China Petroleum Institute,1982,6(2):1-20

[14] 徐世浙.地球物理中的有限元法[M].北京:科技出版社,1994:1-308

Xu S Z.Finite Element method for geophysics[M].Beijing:Science Press,1994:1-308

[15] 王勖成.有限单元法[M].北京:清华大学出版社,2003:472-475

Wang M C.Finite element method[M].Beijing:Tsinghua University Press,2003:472-475

[16] 王月英,宋建国.波场正演模拟中Sarma边界条件的改进[J].石油物探,2007,46(4):359-362

Wang Y Y，Song J G.Improvement of Sarma boundary condition in wavefield forward modeling[J].Geophysical Prospecting for Petroleum,2007,46(4):359-362