弹性波有限差分模拟中自由表面的自适应表达

曹 健,陈景波

(1.中国科学院地质与地球物理研究所,中国科学院油气资源研究重点实验室,北京100029;2.中国科学院地球科学研究院,北京100029;3.中国科学院大学,北京100049)

固体地球表面通常也被称为自由表面,是一个强阻抗突变界面,它的存在使得其附近的地震波波场变得更为复杂,例如会产生强振幅面波(瑞利波和洛夫波),在地表处产生反射纵波和转换横波,以及由于纵横波干涉造成的场地效应等[1]。由于接收到的地震数据大多来自地表或近地表,这些复杂的地震波不可避免地被记录在地震记录中,因此需要在地震波正演数值模拟时考虑自由表面的影响。而自由表面边界条件则是对这一强阻抗突变界面的数学表达。

对于有限元类方法,例如谱元法[2-3]和间断有限元法[4],自由表面边界条件是自然满足的,但这类方法实现起来相对复杂,并且对于复杂模型的计算网格设计较为困难[5]。相比较而言,有限差分法具有算法简单、可灵活用于各种复杂模型的数值模拟、精确和计算高效等优点,目前已被广泛应用于地震波的正演数值模拟。然而,有限差分法中,自由表面的实现需要进行特别处理。为了获得精确并且稳定的自由表面边界实现方式,在过去的几十年中研究人员围绕自由表面边界条件实现方式开展了大量研究。简单地讲,自由表面边界实现方式可分为两大类。第一类是基于弹性波方程的二阶位移表达形式,这种情况下,自由表面边界条件通过位移的空间导数来表示,其形式属于偏微分方程中的Neumann边界条件。因此,可采用不同种类的差分格式来近似自由表面边界处的空间导数,例如,中心差分显式近似格式[6]、单边显式近似格式[7-8]、中心差分隐式近似格式[9]以及复合近似等格式[10-12]。这类自由表面边界处理方式的主要缺陷在于对高泊松比介质进行正演模拟时会产生计算不稳定现象,并且只能达到二阶的计算精度。第二类实现方式是基于弹性波方程的一阶速度-应力表达形式。与第一类方式相比,自由表面边界条件可以直接采用应力为0的形式进行表达,即偏微分方向中的Dirichlet边界条件,这一形式在数学上要相对容易处理一些。其中应用比较广泛的有真空格式[13-14]、应力镜像法[15-17]和一些介质参数修正方法[18-19]。它们的主要优势是对于任意的泊松比值介质都能提供稳定的正演模拟结果,但存在着某些泊松比值情况下计算结果不够准确,需要进行网格加密处理以提高计算精度。

为了解决上述提及的直接针对自由表面边界条件数值近似处理带来的稳定性和精度方面的问题,本文采取了CAO等[20]提出的重新建立自由表面表达方式的想法,从平均介质理论[21]出发,并结合极限思想,从物理角度给出了自由表面的另一种新的自适应数学表达方式,从而更易于在有限差分框架下实现。首先基于平均介质理论,从物理角度给出自由表面处的密度和本构关系应满足的自适应表达方式;然后具体介绍了该自由表面自适应表达方式在速度-应力交错网格有限差分中的实现方式;最后通过一系列数值实验来对比和讨论该自由界面自适应表达方式在不同泊松比情况下的数值模拟效果。

1 自由表面的自适应表达

从物理的角度而言,自由表面是地表以上的空气介质和地下弹性介质之间的一个接触界面,即流体/固体边界。由于地表以上的空气介质,其密度和弹性模量要远远小于地下岩石介质的密度和弹性模量,因此完全可以忽略。根据这一真空近似[22],自由表面还可以被视为真空/固体边界。MOCZO等[21]从平均介质的角度出发,通过弹性模量调和平均和密度参数代数平均的方式来表示各向同性弹性介质内部的介质参数不连续界面,即固体/固体边界。我们将上述平均介质思想引入到自由表面这一强不连续界面的表示上来,从而推导出自由表面的另一种自适应数学表达方式。

假设两个具有不同介质参数的各向同性弹性半空间,其分界面位于z=0处,上层介质的密度和拉梅常数分别为ρ-,λ-和μ-,而下层介质的密度和拉梅常数分别为ρ+,λ+和μ+(图1)。当上层介质为真空时,即ρ-=0,λ-=0和μ-=0,分界面z=0自然就变成了自由表面。根据MOCZO等[21]中提出的介质参数平均方式,分界面z=0处的弹性参数可以表示为:

(1)

式中:κ是体积模量,κ=λ+2μ/3。如果将这一参数平均方式直接应用到自由表面处的介质参数计算上来,即让ρ-=λ-=μ-=0,可以直接得到:

(2)

这里的密度减半可以被理解为自由表面处的质点与内部质点相比更容易发生振动[18],但结果中的拉梅常数为0显然没有意义且不合理。这是因为MOCZO等[21]中最后给出的参数平均方式是针对介质内部的固体/固体不连续界面的,并没有考虑到自由表面的特殊性。因此,我们需要将自由表面的特殊性合理地引入到平均介质的表达方式中。

图1 推导自由表面处本构关系时的弹性参数改变方式

根据MOCZO等[21]的平均介质思想,两个相连的各向同性弹性介质可以用一个平均介质来进行等效,其满足的本构关系为:

(3)

式中:τA和εA分别为平均介质的应力和应变张量,也可以理解为两个介质相连界面处的应力和应变张量,在二维情况下其分量形式为:

(4)

而EA是平均介质的弹性矩阵,同时包含了两种介质各自的拉梅常数,即:

(5)

其中:

(6)

接下来通过逐步改变式(5)中上部介质的物理参数,使得这一固体/固体边界的物理性质逐步近似于自由表面,进而得到可以描述自由表面处应力和应变张量的本构关系。

首先,使上部介质的剪切模量μ=0(如图1),从而上述的固体/固体边界转变为流体/固体边界,其相应的弹性矩阵可以通过将μ-=0代入公式(5)得到,即:

(7)

其中:

(8)

由于空气介质属于流体介质中的一种,所以这一处理方式得到的公式(7)同样可以用来表示空气/固体边界。而自由表面这一概念则是对于地球表面这一空气/固体边界在真空近似下的描述,因此一种比较直接的方式是使公式(7)中上部流体介质的λ=0,来获得真空/固体边界处的弹性矩阵,即自由表面处的本构关系。然而,令λ-=0会使得(8)式中0出现在分母的位置上,进而使得整个弹性矩阵(7)无意义。一种数学上严格且合理的处理方式是利用极限思想,即让λ-趋近于0而不是等于0(如图1),从而可以在极限意义上得到自由表面处的弹性矩阵:

(9)

这里使λ-接近0但不等于0可以很好的表达出地球表面的物理性质,即上部空气介质的弹性模量远远小于地下介质,而使λ-趋近于0则是体现了自由表面这一概念中的真空近似思想。将(9)式代入到(4)式中,我们可以得到自由表面处应力和应变张量满足的本构关系,即:

(10)

此外,引入介质泊松比

(11)

并消去λ,(10)式还可以表示为:

(12)

从中我们可以发现一个有趣的现象。当地下介质是非常致密坚硬的岩石时,即泊松比接近于0,公式(12)在这一极限条件下变为:

(13)

这一结果与MITTET[18]从TI介质近似角度提出的自由表面处理方式的结果相同。而当地下介质为松散的沉积物时,即泊松比接近于0.5时,公式(12)在这一极限条下可以表示为:

(14)

其结果与XU等[19]提出的声波/弹性波边界方法(AEA)相一致。换句话说,MITTET[18]和XU等[19]中的自由表面处理方式分别是新得到的本构关系(公式(10))在小泊松比和大泊松比下的极限情况。这也意味着这一自由表面处的本构关系将适用于更大泊松比范围内的自由表面表达。

因此,我们可以使用(2)式中的密度修正和(10)式中新推导出的本构关系来表达自由表面,并且这一表达方式能够很好地整合到原始的弹性动力学方程中,从而得到自由表面处的弹性波控制方程,即:

(15)

和修正后的本构关系:

(16)

其中ux和uz分别是位移向量沿x方向和z方向的分量。由于(16)式是泊松比自适应的,所以我们将(15)式和(16)式合称为自由表面的自适应表达方式。简单的讲,这一表达方式仅仅是对原始弹性动力学方程中的密度和本构关系的一个修正,因此在计算过程中是不会增加额外的CPU和内存需求。

2 在速度-应力交错网格有限差分中的应用

本节中,我们将给出自由表面自适应表达方式在时间域有限差分数值模拟中的具体实现方式。我们采用显式的交错网格有限差分离散格式,主要是考虑到这种差分格式不会出现高泊松比值时数值模拟精度降低这一问题[20],从而保证自由表面表达方式的泊松比自适应性不会因为数值格式的离散化而受影响。

以VIRIEUX[23]提出的二阶速度-应力交错网格有限差分格式为例,如图2,i,j分别为x方向和z方向格点坐标,自由表面位于j=0处,(15)式和(16)式离散后可得:

(17)

和

(18)

式中:速度分量vx=∂ux/∂t,vz=∂uz/∂t;Δt是时间步长;Δx和Δz分别为x方向和z方向的网格间距。值得注意的是,由于速度垂直分量vz在离散过程中位于自由表面之下的半网格点上(如图2),属于介质内部点,所以(17)式中的第二个式子不需要对密度进行减半修正。此外,(17)式中的应力张量组分需要根据其相对于自由表面的分布情况分为3类进行计算。

图2 自由表面的自适应表达在速度-应力交错网格中的离散方式

第一类,位于自由表面之上(j<0),这种情况下,需要根据真空近似的思想,将应力张量组分全部置0,即:

(19)

第二类,位于自由表面处(j=0),此时应力张量组分需要满足自由表面处的本构关系((10)式),(18)式给出了这一类应力张量组分在时间迭代过程中的计算方式。

第三类,位于自由表面之下(j>0),此时应力张量组分均位于介质内部,不再需要考虑自由表面的影响,可以直接根据这一层介质的本构关系进行计算,即:

(20)

公式(17)、公式(18)、公式(19)、公式(20)只是用来计算自由表面处的速度波场,而模型内部其他位置处的波场仍然采用传统的速度-应力交错网格有限差分格式[23]计算,以这种方式可以很容易地在地震波正演数值模拟中引入自由表面的影响,从而实现全波场的模拟计算。类似地,按照以上3类处理方式将很容易得到这一自由表面表达在高阶有限差分格式中的实现。

3 数值实验与讨论

为了讨论自由表面自适应表达方式在弹性波时间域数值模拟中的计算效果,我们设计了两个数值实验。实验一是模拟均匀半空间中的波场,也就是Lamb问题计算。由于Lamb问题存在解析解[24],因此可以将由自由表面自适应表达方式得出的数值模拟结果与解析解对比,从而评估这一自由表面自适应表达方式的准确性。此外,为了测试其在不同泊松比值情况下的效果,我们设计了4种具有不同泊松比的模型,σ分别为0.200,0.300,0.400和0.489,并将自由表面自适应表达方式计算出的结果与使用MIT-TET[18]和XU等[19]方法得到的结果进行对比。实验二采用了比较接近真实地下结构的Overthrust模型,用来测试这一自由表面自适应表达方式在复杂模型地震波全波场模拟中的可行性。

实验一中介质的纵波速度为5500m/s,密度为1500kg/m3,横波速度分别为3368,2940,2245,807m/s,分别对应于泊松比σ为0.2,0.3,0.4,0.489的情况。模型网格大小为800×200。为了模拟半空间,除顶部的自由表面边界外,其余边界处均使用PML吸收边界条件来消除人工反射。震源是一个单点垂直力,施加在自由表面下的第11个网格层,其震源时间函数是主频为30Hz的雷克子波。在模拟过程中,我们对3种自由表面的处理方式均使用二阶交错网格有限差分格式。此外,为了突出3种数值模拟结果的精度差异,这里不考虑将计算网格加密,仅采用二阶格式体波频散分析得出的采样准则[23],即最小波长内网格点数为10。

图3和图4分别给出了不同泊松比值情况下模拟得到的质点速度水平分量和垂直分量的地震记录,其中红色实线为Lamb问题的解析解结果,蓝色虚线为MITTET[18]方法结果,绿色虚线为XU等[19]提出的AEA方法结果,黑色虚线为自适应自由表面表达方式结果。对比可以看出,3种方法在体波的模拟上都比较精确,而模拟结果的差异体现在瑞利波上。当泊松比较小时,由MITTET方法[18]计算出的地震记录与解析解一致性较好,而XU等[19]方法得出的结果中瑞利波的到达时间较早,与解析解有明显的偏差。随着介质泊松比的增加,XU等[19]方法计算出的地震记录与解析解越来越接近,特别是当σ=0.489时,两者的地震记录基本上完全重合,但使用MITTET[18]方法模拟得到的瑞利波则显示出到时延迟和数值频散现象,且越来越偏离解析解。这意味着这两种方法的模拟效果与泊松比有关,二者的适用性受泊松比值范围的限制。与这两种方法不同的是,自由表面自适应表达方式计算出的结果无论是在小泊松比值还是大泊松比值情况下,都显示出与解析解较好的一致性。为了定量评估3种方法的精度,表1给出了正演模拟地震记录的L2范数误差[14]。可以看出,通过自由表面自适应表达方式计算得到的数值结果,其准确度基本达到了90%,而在小泊松比值和大泊松比值情况下,其结果分别与MITTET[18]和XU等[19]方法得到的结果相近,进而证实了第1节中关于自由表面自适应表达方式在极端泊松比值情况下的讨论结果。

图3 不同泊松比值情况下的Lamb问题合成地震记录(质点速度水平分量vx)对比a σ=0.200; b σ=0.300; c σ=0.400; d σ=0.489

图4 不同泊松比值情况下的Lamb问题合成地震记录(质点速度垂直分量vz)对比a σ=0.200; b σ=0.300; c σ=0.400; d σ=0.489

在实验二中,我们选取Overthrust模型中的部分结构来进行数值模拟测试,这是考虑到在实际情况中,推覆体构造更容易造成地表附近的介质泊松比值变化范围大且不连续。图5展示了计算中使用模型的纵波速度、横波速度、密度和泊松比。在模拟计算中,网格大小为269×1601,网格间距为2.8m。模型顶部边界采用自适应自由表面表达方式进行处理,其余边界使用PML吸收边界条件消除人工反射。与实验一相同,震源为单点垂直力,作用在(x=2990.4m,z=56m)处,其震源时间函数为主频15Hz的雷克子波。图6和图7分别给出了正演模拟得到的不同时刻速度波场的水平分量和垂直分量,从中可以很好的观测到由于自由表面引起的反射纵波和转换横波,以及强振幅瑞利波的产生和传播。此外,通过与纵波速度和横波速度的叠加显示,还可以看出瑞利波在层状介质中的速度频散和泊松比横向突变点处的散射。这些都表明自由表面自适应表达方式的引入能够很好地产生出与自由表面相关且合乎物理规律的波现象,从而证实了其在复杂模型正演模拟中的可行性。

表1 3种自由表面实现方式在不同泊松比情况下的L2误差比较

注:M为Mittet方法结果,X为Xu等提出的AEA方法结果,A为自适应自由表面表达方式结果。

图5 Overthrust模型a 纵波速度; b 横波速度; c 密度; d 泊松比

图6 波场快照(质点速度水平分量vx)a t=0.24s时刻; b t=0.78s时刻; c t=1.05s时刻

图7 波场快照(质点速度垂直分量vz)a t=0.24s时刻,; b t=0.78s时刻; c t=1.05s时刻

4 结论

本文以平均介质理论和极限思想为基础,通过分析自由表面的物理性质,提出了自由表面边界的另一种新的自适应表达方式。简单地讲,这一表达方式是对原始弹性动力学方程在自由表面处的密度和本构关系的修正,因此在数值模拟中不会产生额外的计算量和存储需求。在时间域速度-应力交错网格有限差分中的应用证实了这一表达方式能够很容易地在有限差分数值模拟方法中实现。数值实验和理论分析均表明,这一自由表面自适应表达方式是泊松比自适应的,对于不同大小的泊松比值均能给出较为精确的模拟结果,在空间采样率为最小波长内10个网格点的二阶交错网格差分格式中,其准确度达到90%。此外,推覆体模型的正演模拟证实了这一表达方式在复杂模型正演中的可行性,进而可以应用于实际情况中近地表不均匀风化层内的复杂波场的模拟和分析。