TI介质中P-SV波传播的SBP-SAT模拟

杨在林, 魏夕杰, 孙 铖, 杨 勇

(1.哈尔滨工程大学 航天与建筑工程学院,哈尔滨 150001; 2.哈尔滨工程大学 先进船舶材料与力学工业与信息化部重点实验室,哈尔滨 150001; 3.广西民族大学 建筑工程学院,南宁 530006)

地震模拟方法作为正演过程中的有力工具和策略,也是地震勘探过程中十分重要的一环,它立足于不同介质的参数特性,通过数值计算的方式来模拟波在介质中的传播过程,最终可得到波的传播规律。使用各向异性模型表示具有空间变化的断裂方向和密度、应力分布、复杂层理等的岩石具有特别重要的意义。波动方程数值直接解法以网格化的模型作为基础,将复杂的地形地质条件转化成可以计算求解的形式。它包含了地震波传播过程中的很多信息,因此得到的地震记录往往比较全面真实,模拟精度相对较高。地震波动方程的数值解法主要包括:有限差分法、有限元法、边界元法、谱元法、伪谱法等。

自Alterman等[1]开始使用其进行地震波场模拟至今,有限差分法已发展成为地震数值模拟领域最广泛采用的数值解法之一,具有计算速度快、层状介质模拟能力好、编程简便等优点。有限差分方法也存在需要继续完善改进的地方,比如在模拟波在复杂形状模型边界和内部传播时易产生稳定性问题,时空域长期稳定保证难度较大。这些都是我们期望优化的方向。基于使用能量法证明稳定性的目的,Kreiss等[2]提出了一阶导数的SBP算子基本理论;Carpenter等[3]推导出具有最小模板的二阶SBP算子形式,其内部算子具有与一阶导数内部算子相同的节点数量。Mattsson[4]对SBP算子形式进行修改,使其更加简洁。在处理带边界或界面问题时若使用强施加边界条件,将会破坏SBP算子的性质。得益于Funaro等[5]使用惩戒项的启示,Carpenter等[6]提出了通过SAT方法弱施加边界条件。Mattsson[7]系统地对比了不同方法施加边界条件后对离散系统稳定性的影响,确定采用SAT方法的离散系统更具优势,SBP-SAT离散框架由此确立起来。国内对该方法逐渐重视,Sun等[8-10]在SBP-SAT方法离散框架的基础上搭配完全匹配层作为人工边界,有效压制了边界反射波的产生,获得了稳定的模拟结果。杨在林等[11]推导出了弹性波动方程SBP-SAT离散方程,进一步探究了该方法在波动领域的应用。

经典反射地震勘探通常要求地表水平且地下介质水平层状,但随着勘探领域向更复杂地质区域深入,这些假设条件变得不再适用[12-16]。Igel等[17]将交错网格有限差分法应用于各向异性介质,并且实现了柱坐标系和球坐标系下的数值模拟。不同数值解方法在横向各向同性(transversely isotropy,TI)介质运动学属性研究上各有特色[18-20],对于SBP-SAT方法,在曲线坐标系下对角范数建立的SBP算子能够满足分部求和的性质,进一步的研究由此展开[21]。Carpenter等[22]和Nordström等[23]将SAT方法拓展到曲线坐标下线性问题的多种边界条件及块状界面条件。Virta等[24]建立了适用于复杂形状与非均匀介质的声波SBP-SAT框架,可用于模拟介质不连续模型。孙铖[25]建立了曲线网格矩阵对称型波动方程的多块SBP-SAT方法,并建立了声波SMF体系。SBP-SAT方法与适当的吸收层或边界条件结合并在曲线域中模拟波的传播,公式推导和程序编译得到简化,仿真过程稳定性强,应用前景非常广泛。

本文基于速度-应力形式的弹性波动方程,采用SBP-SAT方法进行了矩阵对称形式的波动方程离散,使用能量法证明其稳定性。对几种TI介质模型进行数值模拟,得到波场信息来分析其中P-SV波的传播规律,并对比COMSOL结果验证SBP-SAT方法的正确性。

1 波动方程SBP-SAT方法构建

1.1 弹性波动方程的矩阵对称化

常见TI介质模型示意图,如图1所示。倾斜横向各向同性介质速度-应力形式的P-SV波波动方程(体力为零)为

图1 常见TI介质模型示意图

(1)

式中:ρ为弹性体密度;σij为应力;vx和vy分别为质点速度在x轴和y轴上的分量。为表述方便仅考虑绕z轴旋转情形,设极化角为θ,刚度矩阵元素dij为

d11=(C11cos2θ+C13sin2θ)cos2θ+

(C13cos2θ+C33sin2θ)sin2θ+C44sin2(2θ)

d13=(C11cos2θ+C13sin2θ)sin2θ+

(C13cos2θ+C33sin2θ)cos2θ-C44sin2(2θ)

d15=0.5(C11cos2θ+C13sin2θ)sin(2θ)-

0.5(C13cos2θ+C33sin2θ)sin(2θ)-

C44sin(2θ)(cos2θ-sin2θ)

d33=(C11sin2θ+C13cos2θ)sin2θ+

(C13sin2θ+C33cos2θ)cos2θ+C44sin2(2θ)

d35=0.5(C11sin2θ+C13cos2θ)sin(2θ)-

0.5(C13sin2θ+C33cos2θ)sin(2θ)+

C44sin(2θ)(cos2θ-sin2θ)

d55=0.25[C11sin(2θ)-C13sin(2θ)]sin(2θ)-

0.25[C13sin(2θ)-C33sin2θ)]sin(2θ)+

C44(cos2θ-sin2θ)2

(2)

将波动方程表示为矩阵形式

(3)

(4)

M=PTΛ-1/2P,M-1=PTΛ1/2P

(5)

代入(3)中得到

(6)

进一步,得

(7)

至此我们得到了SMF的弹性波动方程

(8)

传统有限差分法在模拟曲线域时大量角点会影响模拟结果的稳定性。将曲线坐标系映射到直角坐标系并划定合适的网格形式能更好模拟层状地质结构。

将曲线域Ω∈(x,y)转换成直角域Ω′∈(ξ,η),如图2所示。坐标映射关系为

图2 坐标变换

[x(ξ,η),y(ξ,η)]↔[ξ(x,y),η(x,y)]

(9)

为实现这一映射关系需要使用雅可比矩阵行列式和链式法则

(10)

曲线域的参数向量方向导数为

JUx=Uξyη+UηyξJUy=Uξxη+Uηxξ

(11)

利用式(10)可以得到弹性波动方程式(8)在曲线域的矩阵对称形式为

Ut=AξUξ+AηUη

(12)

由于Ax和Ay为对称矩阵,所以Aξ和Aη也为对称矩阵。参数矩阵为

(13)

拓展SBP算子得

(14)

Iξ和Iη分别为(Nξ+1)阶矩阵和(Nη+1)阶矩阵,拓展系数矩阵得

(15)

边界向量为

EW=e1ξ⊗Iη⊗I5EE=erξ⊗Iη⊗I5ES=Iξ⊗elη⊗I5EN=I1ξ⊗erη⊗I5

(16)

求解域边界定义为W,E,S,N,得到SBP-SAT方法建立的矩阵对称型弹性波动方程的离散形式为

Ut=AξDξU+AηDηU+SATW+SATE+SATS-SATN

(17)

式中,SATW,SATE,SATS,SATN的表达式分别为

(18)

1.2 适定的边界条件与半离散近似

通过特征边界条件的推导能够得到传入和传出求解域的物理量。波动方程变形得到

(19)

(20)

对于SMF框架来说,设置人工边界的各种边界性质只需要选用相应的反射系数即可。

p=Rq

(21)

反射系数的定义域为[-1,1],取值为0时即无反射边界条件,取值为1时即自由边界条件。传入和传出的特征变量关系可以定义不同的边界条件,也意味着更大的适用范围。边界条件的通用框架建立起来。

1.3 系统能量分析

结合1.2节中得到的边界条件式(20),对弹性波动方程在连续域进行能量分析。左乘UT后与其转置相加,并在求解域上进行积分,得到

(22)

式(22)等价于

(23)

边界项BTW,BTE,BTS,BTN分别为

BTW=-∬W[UT(Ax)U]dSBTE=∬E[UT(Ax)U]dSBTS=-∬S[UT(Ay)U]dSBTN=∬N[UT(Ay)U]dS

(24)

边界条件为无反射边界条件时,BTW,BTE,BTS,BTN均为非正值,说明系统的能量有界,证明矩阵对称型弹性波动方程框架满足稳定性的要求。弹性波动方程的矩阵对称形式使能量分析过程简化,并建立起了通用框架。

离散化的弹性波动方程的证明过程与上述过程类似,当满足边界条件式(20)时,式(18)恒等于0。

将式(17)左乘UTHξη再与该式子转置相加后,得

(25)

其中,

边界项值均为非正,说明矩阵对称型弹性波动方程的离散框架满足稳定性要求。

2 算例分析

2.1 倾斜横向各向同性介质模型SBP-SAT方法模拟

倾斜横向各向同性(tilted transverse isotropic ,TTI)介质是具有倾斜角度对称轴的TI介质形式,其刚度矩阵可对垂直横向各向同性(transverse isotropy with a vertical axis of symmetry,VTI)介质弹性系数矩阵进行Bond变换而得到,如图3所示。

图3 多块模型示意图

建立的模型尺寸1 980 m×1 980 m,有3×3个子计算域,每个计算域的网格数为200×200,采用4阶精度传统SBP算子,时间步长Δt取1 ms。定义模型的弹性系数分别为C11=28.00 GPa,C13=8.00 GPa,C33=15.00 GPa,C44=3.00 GPa,极化角取45°。初始条件设置为

g(x,y)=exp{-α[(x-x0)2+(y-y0)2]}

(26)

α影响波前面,取0.002目的是使波场效果更清晰明显。单炮记录位置选择在深度825 m处。

获得的结果包括波场快照和单炮记录,在图4的波场中存在波和波。同时波场没有明显的数值频散现象,展现出SBP-SAT方法数值模拟的保真性和精准性。

图4 极化角45°时的单炮记录和x,y方向速度幅值

在TTI介质中由于波的偏振方向与传播方向存在夹角,除特定传播方向外不存在纯横波或纯纵波,需要使用克里斯托弗方程求解偏振向量,相关波场分离可以作为下一步的研究方向。分析图4可知,介质中波的传播速度受到介质参数影响不能保持相同,波场快照中外侧的为波,内侧的为波,并且波的波前面相较于算例2.2对应位置更圆,波波前面角度变化并出现波面尖角现象。横纵波场是耦合在一起的。

2.2 分层VTI介质模型SBP-SAT方法模拟

模型的尺寸、子计算域数量和网格数与2.1节中相同,采用6阶精度传统SBP算子,时间步长取1 ms。多块模型分上下两层,上层0～1 320 m,对应子计算域2、3、5、6、8、9;下层1 320～1 980 m,对应子计算域1、4、7。极化角均取0°,构成分层的VTI介质模型。

初始条件和单炮记录位置均不变。

表1 双层VTI介质模型参数表格

本算例主要为分析在分层界面处P-SV波的转换和传播规律。观察图5可以看到在同一介质内部,波前面清晰连贯,但在介质分界处有明显的波前面断裂和速度差异,且上下两部分形成明显完整的反射波、透射波和转换波。

图5 单炮记录和x,y方向速度幅值

为了震相的简单标注,在此规定:或为波或波;角标1、2分别为上、下两层;上标为上行波;下标为下行波。分析图5可以得到数字:1为直达qP波;2为直达qSV波;3qP1qp1为反射波;4qP1qS1为反射转换波;5qSV1qSV1为反射波;6qSV1qSV2为透射波;7qP1qSV2为透射转换波;8qP1qP2为透射波。

从图5可知,随时间增大,波前面也在增大,波和波耦合在一起传播,只在某些特殊方向上独立传播。当波遇到波阻抗分界面时会发生转换,波不仅会产生反射波和透射波,还会产生反射转换波和透射转换波。波也有同样的现象,其波场变化比较明显,由于介质的速度各向异性导致波的速度随传播方向变化而变化,产生了波面尖角。

另外波前面到达边界后没有明显的反射,证明无反射边界的效果良好。

2.3 含有裂缝的VTI介质模型SBP-SAT方法模拟

不同地质年代沉积形成的地层可近似为TI介质,而其中裂缝、断层是普遍存在的,使用数值模拟方法对含裂缝的VTI介质开展地震动模拟研究具有丰富的工程背景。

建立的模型尺寸、子计算域数量和每个计算域的网格数均不变,采用6阶精度传统SBP算子,时间步长Δt取1 ms。裂缝位于子计算域4～5,裂缝曲线设置为正弦函数,子计算域5为曲线域,如图6所示。定义弹性系数取C11=28 GPa,C13=8 GPa,C33=15 GPa,C44=3 GPa,极化角取0°。

图6 多块模型示意图

初始条件施加在方向上,其他条件不变,如图7、图8所示。

图7 x方向速度幅值

图8 y方向速度幅值

SAT弱施加的特点使得边界设置富有多样性,在裂缝、空腔等类型边界条件建立上,不必采用介质填充等相对繁琐的方式,而是可以通过设置反射参数这种简单的方式来完成。本算例中子计算域4和5边界的反射系数均设置为1,即设为自由边界。观察图7,曲线域中波的传播连续清晰,波场效果稳定,同时裂缝边界的反射作用也很明显。

进一步,对相邻两子计算域网格、节点位置不对应的情形,可以采取子计算域各自建立更优的节点分布的策略,这大大拓宽了方法的应用范围。

3 准确性验证和计算效率

为求简明,建立均匀各向同性介质模型,尺寸设置为3 000 m×3 000 m,横波波速为2 000 m/s,纵波波速为1 100 m/s,密度为1 100 kg/m3。初始条件为

(27)

采用4阶精度传统SBP算子,时间步长为0.001 s,子计算域网格划分为200×200。通过节点(1 500 m,1 000 m)的x方向求解速度幅值,与COMSOL软件使用间断伽辽金法计算的结果进行比较,观察幅值差随时间变化。

如图9和图10所示,结果吻合度较好,能够验证结果的准确性。

图9 x方向速度幅值

图10 x方向的节点幅值差

采用三种网格大小来计算CPU时间,模型采用2.1节的模型,网格数为150×150,300×300,600×600,迭代500次,分别得到的CPU计算时间为6.564 s,21.368 s,110.048 s。

4 结 论

本文使用SBP-SAT方法,对不同类型的横向各向同性介质模型进行了波场数值模拟研究。结果表明,SBP-SAT方法表现出了很好的适用性。在模拟分层介质和含裂缝的情形时,既能够保证计算精度又能够保证稳定性,进一步优化子计算域划分方式下有继续提高计算效率的潜力。SBP-SAT方法同样适用于其他介质中波动方程的求解,推导过程与本文讨论的方法相近。