时间域广义2.5D一阶波动方程曲网格有限差分法数值模拟

杨尚倍 白超英 ZHOU Bing

(①中国地质调查局西安地质调查中心，陕西西安 710054；②长安大学地质工程与测绘学院地球物理系，陕西西安 710054；③哈利法科学与技术大学地球科学系，阿布扎比 2533)

0 引言

地震波场数值模拟不仅是了解地震波传播规律和解释地震数据的重要工具，也是地震成像的基础。随着计算机技术快速发展，各种地震波场数值模拟算法日趋完善，如有限差分法[1]、有限元法[2]、边界元法[3]、伪谱法[4-5]、谱元法[6]、有限体积法[7]等。尽管3D波场数值解相对于2D更符合实际观测记录，然而巨大的内存需求和运算量使得实际3D大尺度波场数值模拟及波形反演面临计算资源的巨大挑战。因此，目前更多的实际应用则在2D或小尺度的3D空间开展，如逆时偏移成像[8]和全波形反演[9]。

虽然2D数值模拟相较于3D高效，但引入了线源假设代替实际的点源激发。这种线源假设意味着震源激发大小沿走向不变，模拟的波场也是2D，其结果很难与实际观测的3D地震记录相匹配[10-11]。当然可采用3D到2D地震数据转换方法，将实际的点源数据(3D)转换成线源数据(2D)，但目前这些转换方法只适用于简单的地质模型，例如均匀介质和水平层状介质，并且需要满足远场近似和无波场重叠等苛刻条件[12]，不利于实际复杂介质的成像。

鉴于此，为了克服2D波场模拟中线源假设带来的失真和3D波场模拟中效率低下的不足，人们提出了2.5D波场模拟方法[13]。2.5D地震波场数值模拟方法假设地质模型沿走向方向模型属性不变，因此可以沿走向方向做傅里叶变换将3D问题简化为重复的2D问题，从而大大提高了计算效率，且与3D解析解高度吻合[14]。如：Song等[13]在频率域通过有限差分法重构了3D地震波场；Zhou等[15]给出了精确的2.5D吸收边界条件，并采用有限元法得到了频率域2.5D波场数值解；Maokun等[16]则利用近似最佳正交有限差分算法进行了2.5D弹性波数值模拟；Takenaka等[17]给出了时间域平面波入射下2.5D弹性波动力学方程的求解过程；Novais等[18]改进了时间域2.5D地震波场有限差分数值模拟方法的波数采样方法。Xiong等[19]首先将空间域弹性波方程转换为波数域弹性波方程，然后采用交错网格有限差分求解波数域方程。

然而，上述时间域或频率域2.5D地震波场数值模拟方法只针对单一的弹性、黏弹性或各向异性介质，并没有一种较为普适的算法，即可同时满足不同介质属性(弹性、黏弹性、各向异性)和不同边界条件(声波自由地表，固体自由地表和固—液边界)下的2.5D波场数值模拟算法。本文提出一种适合复杂多种介质混合模型下时间域广义2.5D一阶波动方程及数值解法，适用于混合(声波、弹性各向同性、弹性各向异性)介质和各种边界条件(声波自由地表、固体自由地表和固—液边界)的地震波场数值模拟。

另一方面，为了适应复杂地表起伏模型中的数值模拟，本文采用曲网格进行模型剖分，并采用低频散、低耗散的同位网格MacCormack有限差分法格式[20]离散时间域广义2.5D一阶波动方程，从而得到高精度的数值解。最后，应用数值模拟算例验证了本文2.5D数值模拟方法的正确性和有效性。

1 时间域2.5D广义一阶速度—应力 波动方程

首先建立广义2.5D一阶速度—应力波动方程，并分析该方程的普适性问题。2.5D地震波场数值模拟是3D地震波场数值模拟的一种特殊情况，其假设地质模型沿走向(y方向)具有不变的物理性质，则2.5D地震波动方程可以通过对3D地震波动方程沿y方向做傅里叶变换得到[21]。

1.1 弹性各向异性介质波动方程

三维弹性各向异性介质中，地震波的控制方程[22]为

(1)

σ=cu

(2)

(3)

和

(4)

(5)

式中上标“(S)”表示弹性各向同性或弹性各向异性的固体介质。式(3)和式(4)可简写为矩阵形式

(6)

上式为任意弹性各向异性介质中的2.5D地震波控制方程，系数矩阵A(S)、B(S)、C(S)的具体表达式见附录A。

1.2 声波介质波动方程

在声波介质中，应力满足切向应力为零和法向应力的轴向分量等于液体压力

σij=-δijP

(7)

P可表示为

P=-Kuj,j+s(t)δ(x-xs)

(8)

式中：K为体积模量；s(t)为xs处的震源函数。将式(7)代入式(1)，可得

(9)

对式(8)求时间导数，可得

(10)

根据式(9)，式(10)可以表示为

(11)

对式(9)和式(11)做关于y的傅里叶变换，可得

(12)

(13)

定义

(14)

可得式(12)和式(13)的矩阵形式

(15)

式中：上标“(W)”表示声波介质；系数矩阵A(W)、B(W)、C(W)的具体表达式见附录B。上式为声波介质中的2.5D地震波控制方程。

1.3 声波自由地表

声波自由地表边界条件需要满足地表处液体压力P=0，代入式(12)和式(13)，可得

(16)

(17)

式中：上标“(AW)”表示声波自由地表；系数矩阵A(AW)、B(AW)、C(AW)分别为

(18)

(19)

C(AW)=0

(20)

对比式(18)～式(20)与附录B，可以发现式(17)是式(15)的特例。

1.4 固体自由地表

固体自由地表边界条件需要满足地表处法向应力为0，即

σ·n=0

(21)

展开为

(22)

式中n=(n1，n2，n3)表示界面的法向向量。将式(22)写成矩阵形式

(23)

将上式写为子矩阵格式

Γ2σ2=-Γ1σ1

(24)

式中

(25)

(26)

定义

(27)

式(23)可改写为

σ2=Sσ1

(28)

由上式可知，应力分量σ2可由σ1和表面法向分量n得到，只有应力分量σ1未知。假设自由地表由z=z0(x)给出，其坡度为tanα=z′0(x)，法向分量为n=(-sinα,0,cosα)，则式(27)变为

(29)

根据式(28)和式(29)可得

(30)

将上式代入式(3)，可得

(31)

由式(4)可得

(32)

定义

(33)

则式(30)和式(31)可简写为

(34)

1.5 固—液边界

在固—液边界上，例如由z=z0(x)定义的二维海底界面，边界条件需要满足法向应力和法向速度分量连续，并且切向应力分量为0，即

(35)

(36)

(37)

(38)

将式(38)代入(37)，可得

(39)

(40)

(41)

为了表示简单，假设固液边界上没有施加震源，将式(4)和式(13)代入(41)，可得

(42)

上式对于任意ky成立，可得

(43)

由式(35)可知

(44)

将式(42)代入式(4)，可得

(45)

将式(39)和式(40)代入式(3)，可得

(46)

式中ρS为固体密度。由式(12)可知

(47)

式中ρW为液体密度。定义

(48)

将式(45)～式(47)写成矩阵形式

(49)

至此，获得了2.5D广义一阶速度—应力波动方程

(50)

式中的各个参量可根据不同介质或不同边界条件而有所不同。

2 曲线网格有限差分算法

2.1 曲线坐标系与直角坐标系的转换关系

为了使网格线与物理空间的边界线重合，避免传统矩形网格由于阶梯状近似产生的虚假散射问题，应用曲线网格(贴体网格)对地质模型进行剖分，然后将实际的非规则物理空间变换到一个规则的计算空间[23]。假设物理空间的坐标(x,z)和计算空间的坐标(ξ,η)之间映射关系(图1)为

图1 曲线网格和计算网格的变换关系示意图

(51)

由链式求导法则可得

(52)

整理可得到坐标转换系数

(53)

式中J=x,ξz,η-x,ηz,ξ。为了避免网格非均匀引发的虚假震源项，x,ξ、x,η、z,ξ、z,η均采用数值方法计算，所采用的数值法与下文求解波动方程空间导数的数值算法相同。

2.2 曲线坐标系中的广义一阶速度—应力波动方程

根据求导的链式法则，广义一阶速度—应力弹性波方程的矩阵形式在曲线坐标系中的表达式为

(54)

Hixon等[24]采用Tam等[25]提出的DRP (Dispersion Relation Preserving)方法对同位MacCormack格式的系数进行了优化，得到了低频散、低耗散的DRP/opt MacCormack 格式，本文采用该格式有限差分算子离散方程。

同位网格MacCormack格式的两个偏心差分格式——向前和向后差分格式分别为

(55)

(56)

差分系数为

(57)

采用DRP/opt MacCormack有限差分格式，式(54)空间导数可以离散为

(58)

在时间方向上采用四阶Runge-Kutta法[26]，结合DRP/opt MacCormack有限差分算子，波场可以从时间步长第N步更新到第N+4步[20]，具体过程如下。

(59)

(60)

(61)

(62)

式中：α1=0.5；α2=0.5；α3=1.0；β1=1/6；β2=1/3；β3=1/3；β4=1/6。

3 震源、人工边界和波数域采样

3.1 震源和吸收边界

本文在网格点上施加集中力源，并使用Ricker子波做为震源

(63)

式中：震源中心频率fc=30Hz；延迟时间t0=30ms。为了消除人工边界产生的边界反射，本文采用新近提出的广义刚度衰减法(The Generalized Stiffness Reduction Method，GSRM)[27]，吸收效果优于最佳匹配层方法，适用于各种介质，如声波介质、各向同性和各向异性介质。

3.2 波数域采样方法

(64)

(65)

式中vij表示施加i方向集中力源时得到的j方向速度分量。

(66)

4 2.5D与3D数值模拟的精度和效 率分析

表1 介质弹性参数

图2为三个地表水平均匀介质模型的2.5D数值解的波场快照，从图中可以清晰地看到波场的P、S波波前。为了验证本文方法模拟得到的2.5D数值解的计算精度和正确性，将2.5D数值解与对应的3D解析解及3D数值解进行对比(图3)，三者匹配良好，验证了本文算法的计算精度和正确性。

图2 地表水平均匀介质模型2.5D模拟的波场快照(a)声波介质的P分量，0.15s；(b)、(c)分别为各向同性介质的vx和vz分量，0.09s；(d)、(e)分别为VTI-1介质的中vx和vz分量，0.08s

图3 地表水平均匀介质模型2.5D数值解与3D数值解、3D解析解的对比(a)声波介质中P分量；(b)、(c)分别为各向同性介质的vx和vz分量；(d)、(e)分别为VTI-1介质的vx和vz分量

2.5D与3D数值模拟计算机内存占用和CPU运行时间的对比如表2所示，可以看出，2.5D数值模拟所需CPU时间和计算机内存都远小于3D，特别是计算机内存。所以本文的2.5D波场数值模拟方法则具有较高的实际应用价值，可以直接应用于地震数据反演方法。

表2 2.5D与3D数值模拟用时和内存占用对比

5 复杂模型的波场数值模拟

5.1 起伏地表均匀模型

应用起伏地表声波、弹性各向同性和弹性VTI-1均匀介质模型(图4)验证本文施加自由地表边界算法的正确性。网格剖分间距Δx=Δz=2m，震源位于(500m，400m)，检波器位于(300m，500m)，模型参数见表1。

图4 模型网格剖分示意图红星和蓝三角分别为震源和检波器

图5、图6、图7分别为起伏地表声波、弹性各向同性和弹性VTI-1均匀介质模型的2.5D数值解的vx和vz分量波场快照，从中可以清楚地看到直达P波及自由地表产生的反射波。图8、图9、图10为三个起伏地表模型2.5D与2D数值解的波形对比，其中2D数值解是令式(50)中ky=0后采用本文相同方法获得。2.5D数值解和2D数值解之间不仅存在着明显的振幅差异，还存在相移。

图5 起伏地表声波均匀介质模型2.5D数值解vx(a)和vz(b)分量在0.4s的波场快照

图6 起伏地表各向同性均匀弹性介质模型2.5D数值解vx(a)和vz(b)分量在0.32s的波场快照

图7 起伏地表均匀VTI-1介质模型2.5D数值解vx(a)和vz(b)分量在0.32s的波场快照

图8 起伏地表声波均匀介质模型2.5D与2D数值解对比(a)vx分量；(b)vz分量；(c)图a的归一化显示；(d)图b的归一化显示。下图为上图方框的放大显示

图9 起伏地表均匀各向同性弹性介质模型2.5D与2D数值解的对比(a)vx分量；(b)vz分量；(c)图a的归一化显示；(d)图b的归一化显示。下图为上图方框的放大显示

图10 起伏地表均匀VTI-1介质模型2.5D与2D数值解的对比(a)vx分量；(b)vz分量；(c)图a的归一化显示；(d)图b的归一化显示。下图为上图方框的放大显示

5.2 固—液耦合模型

建立一个正弦起伏固—液耦合界面模型(图11)验证本文的2.5D数值模拟方法适用性。模型尺寸为1000m×1000m，网格间距为2m，第一层为声波介质，第二次为各向同性介质(参数见表1)。震源位于(500m，300m)，5个检波器等间距布置于(300m，0m)与(300m，400m)之间。图11为2.5D数值解的0.28s时刻波场照，从图中可以观察到直达P波经过固—液界面产生了清晰的透射P波和S波；图12是2.5D与2D数值解的波形对比，可以看出2.5D与2D数值解之间同样存在着明显的振幅差异和相移现象。

图11 固—液耦合模型2.5D数值解vx(a)和vz(b)分量在0.28s的波场快照图中起伏曲线为正弦起伏界面

图12 固—液耦合模型2.5D与2D数值解波形对比(a)vx分量；(b)vz分量

5.3 三层混合模型

为了检验本文2.5D波场数值模拟方法对于复杂介质的适用性，选择一个三层模型(图13)进行实验。该模型第一层为声波介质(水)，第二层为弹性各向同性介质，第三层为VTI弹性介质(VTI-2)，各层弹性参数见表1。震源位于(1000m，450m)；三个检波器分别位于(500m，200m)、(500m，500m)、(500m，800m)，每层各有一个。图14为三层模型2.5D数值解vx分量在不同时刻的波场快照，可以看出清晰的反射波、转换波、透射波和多次波。由图14a、图14b可以看出，当P波和S波到达固—液边界时，透射和转换后的P波出现在第一层(声波介质)中，而反射和转换后的P波和SV波在第二层(各向同性介质)中向下传播；在第三层(VTI-2介质)中出现P波和S波穿过固—固界面产生的透射和转换波。由图14c～图14f可以发现，当波传播到自由表面产生了反射的P波，在第二层介质中出现了清晰的多次波。图15、图16、图17分别为三个检波器的2.5D与2D模拟记录对比，可以看出，无论检波器位于何种介质中，二者存在明显的振幅差异和相移。

图13 三层模型示意图星和三角分别表示震源和接收器点位置

图14 三层模型2.5D数值解vx分量在不同时刻的波场快照(a)0.12s；(b)0.20s；(c)0.28s；(d)0.36s；(e)0.44s；(f)0.52s

图15 三层模型第一个检波器的2.5D与2D数值模拟记录对比(a)vx分量；(b)vz分量；(c)图a的归一化显示；(d)图b的归一化显示。下图为上图方框的放大显示

图16 三层模型第二个检波器的2.5D与2D数值模拟记录对比(a)vx分量；(b)vz分量；(c)图a的归一化显示；(d)图b的归一化显示。下图为上图方框的放大显示

图17 三层模型第三个检波器的2.5D与2D数值模拟记录对比(a)vx分量；(b)vz分量；(c)图a的归一化显示；(d)图b的归一化显示。下图为上图方框的放大显示

三层混合模型模拟结果表明，本文方法适用于复杂介质2.5D地震波数值模拟。

6 结论

本文通过推导、整理得到了一个矩阵形式的2.5D广义一阶速度—应力波动方程，该方程不仅适用于声波介质和一般各向异性介质，而且适用于不同的边界条件(声波自由地表、固体自由地表和固—液边界条件)。通过定义不同的波场矢量，并根据模型离散点介质参数的变化给出了三个系数矩阵，在压力或应力矢量上施加点源后采用数值方法求解，可以用一个计算机程序实现复杂介质中的2D或2.5D地震波场数值模拟。为了更好地贴合起伏地形，本文选择曲线网格有限差分法求解波动方程，数值模拟实验表明，在不同的介质中本文方法模拟得到的2.5D数值解与3D解析解非常匹配，并且适用于不同的边界条件。除此之外，本文验证了2.5D数值方法的计算效率远优于3D数值方法。这种2.5D数值模拟方法可以直接用于实际点源地震数据的逆时偏移成像。

附录A 弹性各向异性介质中波动方 程的系数矩阵

将式(3)、式(4)写成矩阵形式，则式(6)中的系数矩阵可表示为

(A-1)

(A-2)

(A-3)

式中

(A-4)

(A-5)

(A-6)

(A-7)

(A-8)

(A-9)

(A-10)

(A-11)

(A-12)

(A-13)

(A-14)

附录B 声波介质中波动方程的系数 矩阵

将式(12)和式(13)写成矩阵形式，则式(15)的系数矩阵为

(B-1)

(B-2)

(B-3)

附录C 固体自由地表中波动方程的 系数矩阵

将式(31)和式(32)写成矩阵形式，则式(34)系数矩阵可表示为

(C-1)

(C-2)

(C-3)

式中：

(C-4)

(C-5)

(C-6)

(C-7)

(C-8)

(C-9)

附录D 固—液边界中波动方程的系 数矩阵

将式(45)～式(47)联合写成一个矩阵形式，可导出式(49)的系数矩阵，为

(D-1)

(D-2)

(D-3)

式中：

(D-4)

(D-5)

(D-6)

(D-7)

(D-8)

(D-9)