各向异性介质中弹性波方程的有限差分格式及其稳定性条件

(石油大学地球科学系，北京100083)

腾张基温中介

(中国科学院地球物理研究所，北京100101)

在勘探地震学中，快速有效地模拟地震波在各向异性介质中的传播具有重要意义。算法的稳定性分析对于加快计算速度是非常必要的。本文首先用矩阵和向量来描述波动方程，针对二维和三维一般各向异性介质中的弹性波方程，提出了一种快速且消耗内存较少的有限差分方法。然后，系统地研究了二维均匀和非均匀各向异性波动方程有限差分格式的稳定性条件。此外，还给出了有限差分法在某些特殊各向异性条件下的稳定性公式。最后，本文还研究了三维有限差分格式的稳定性。

弹性波方程；有限差分；稳定各向异性介质

1简介

地震波传播的数值模拟在地球科学中具有重要意义。在各向异性地震模拟的各种方法中，基于Kennett的研究工作。Kosloff等人[13]用傅立叶方法模拟地震各向异性。而陈[7]采用有限元法模拟非均匀各向异性。虽然有限差分法已被广泛用于模拟各向同性介质中的弹性波，但用这种方法模拟地震各向异性并不常见。Tsingas等人[16]开发了一种利用有限差分算子求解横观各向同性介质中偏微分方程的模拟算法，该算法基于MacCormack型分离格式[2]。Faria等人[8]基于交错网格格式[17]，采用有限差分算法模拟二维横观各向同性问题。最近，Igel等人[9]给出了一种基于卷积算法模拟一般地震各向异性的有限差分算法。

速度快、存储量少是有限差分法的优点。随着大规模波场模拟的需要和大规模并行计算的发展，复杂介质或高维(二维和三维)模型的有限差分地震模拟不可避免。虚拟谱方法是非常流行的，因为它的空间算子精确地达到奈奎斯特频率。但是虚拟谱法需要傅里叶变换，所以模拟三维各向异性非常耗时。同时，采用傅立叶变换意味着每个网格点与其他点相互作用。从某种意义上说，这与动力学的局部弹性性质不一致。因此，在设计地震模拟的有限差分格式时，有必要考虑差分算子的局部性。另一方面，考虑差分算子的局部性对提高算法的并行性非常重要。由于最近的网格点之间的信息交换最快，模拟各向异性大尺度模型的波场是可行的。

基于上述原因，本文提出了一种快速的、占用内存较少的各向异性介质中地震波传播的有限差分方法。其实这是算法[10，18]的推广。

地震波传播的数值模拟通常采用时间推进算法。为了保证算法的稳定性，时间增量受到算法稳定性条件的限制。选择合适的时间步长既能保证数值计算的稳定性，又能加快计算速度。否则，不仅会产生非物理的数值振荡，甚至会导致错误的结果。合理时间增量的选择取决于差分格式和描述介质特性的介质参数或速度，换句话说，取决于差分格式的稳定性条件。因此，有限差分格式的稳定性在数值计算中至关重要。虽然在文献[10，18]中对这一问题进行了一些特殊情况的研究，但他们并没有对一般各向异性问题的有限差分格式及其稳定性进行详细系统的研究。我们的目的是给出该问题的一般有限差分格式及其稳定性条件，并进一步给出一些特殊各向异性条件下的稳定性公式。结果表明，在各向同性条件下，我们的结果与Aboudi[1]的结果一致。

2有限差分格式

2.1三维各向异性

各向异性介质中的运动平衡方程可以是以下形式:

岩石圈结构和深部作用

其中:ρ为介质密度；Fx、fy、fz分别代表力源在x、y、z方向的分量；Ux、uy和uz分别表示x、y和z方向的位移分量。

应力应变关系为，其中弹性参数矩阵，和ci，j=cj，I，I，j=1，2，…，6；

应力矩阵σ = (σ xx，σyy，σzz，σyz，σxz，σxy)t；

应变矩阵ε = (ε xx，εyy，εzz，εyz，εxz，εxy)t；

应力和位移之间的关系是:

岩石圈结构和深部作用

制造

岩石圈结构和深部作用

岩石圈结构和深部作用

那么等式(1)可以写成:

岩石圈结构和深部作用

显然，A，E，Q，B+D，C+G，F+H是实对称矩阵。在不考虑源项的情况下，使用有限差分近似方程(2)可以得到以下有限差分格式:

岩石圈结构和深部作用

其中:δx、δy、δZ分别代表x、y、Z方向的空间步长，δt代表时间步长，

岩石圈结构和深部作用

2.2二维各向异性

类似地，二维各向异性介质中的波动方程可以写成:

岩石圈结构和深部作用

并且有以下不同的格式:

岩石圈结构和深部作用

3稳定性条件

3.1均匀介质中二维有限差分格式的稳定性条件

格式(5)在均匀介质的情况下可以简化为:

岩石圈结构和深部作用

根据Richtmyer和Morton的稳定性理论分析，我们作出

岩石圈结构和深部作用

其中U0是常数向量，等式(6)变成:]]

2I—2Pλ+I)U0=0

其中I表示单位矩阵，

岩石圈结构和深部作用

其中:，，，

如果系数行列式为0，则:

岩石圈结构和深部作用

定理1差分格式(6)在下列条件下是稳定的:

岩石圈结构和深部作用

其中函数F(α，β)和α，β为:

岩石圈结构和深部作用

岩石圈结构和深部作用

其中k = δ z/δ x

根据差分格式(6)的稳定性证明方程(7)中的λ满足‖λ‖≤1。

根据(7)和引理2(见附录)，有以下不等式:

岩石圈结构和深部作用

从A，Q和C+G的对称性可知，矩阵也是对称的，所以由引理3(见附录)和不等式(8)，有

岩石圈结构和深部作用

从矩阵A，Q和C+G的元素的非负性，我们可以使0≤α，要使(9)成立，我们只需要

岩石圈结构和深部作用

制造

岩石圈结构和深部作用

求偏导数

岩石圈结构和深部作用

根据波动方程的特点，有

‖A‖0，Q‖0，C+G‖0

所以(0，0)和(，)可能是函数的极值点。很明显，

f(0，0)=0

岩石圈结构和深部作用

我们来讨论(α，β)≦(0，0)或(，)的情况:

通过(11)和(12)，如果，那么

岩石圈结构和深部作用

其中α和β可能是f(α，β)的极值点，所以稳定性条件为:

岩石圈结构和深部作用

否则就是函数的最大值点，稳定性条件为:

岩石圈结构和深部作用

特别地，当δ x ~ δ z时，有以下简化的稳定性条件:

岩石圈结构和深部作用

其中α和β由(13)和(14)确定，k=1。

3.2非均匀条件下的稳定性条件

在不均匀的情况下，很难确定差分格式(5)的稳定性条件。但根据偏微分方程的数值方法理论[15]，我们可以用“冻结系数”法[14]来分析其稳定性。此外，我们给出了非均匀介质情况下差分格式(5)的稳定性条件。

实际上，如果介质的参数函数是连续有界的，我们一般可以把一个小的计算区域看成是一致的，然后差分格式(5)可以退化成格式(6)，这样就可以得到小区域内的稳定条件也是一致的。此外，根据介质参数函数的连续有界性，我们可以获得差分格式(5)的稳定性条件如下:

如果是这样的话

岩石圈结构和深部作用

否则，

岩石圈结构和深部作用

其中H(α，β)，α和β为:

岩石圈结构和深部作用

3.3差分格式在某些特殊介质中的稳定性条件

显然，通过计算格式(6)中矩阵A、Q和C+G的范数，可以分别获得各向同性和横向各向同性情况下的稳定性条件如下:

均质性

岩石圈结构和深部作用

或者

岩石圈结构和深部作用

其中λ和μ是拉梅常数，是P波速度。

显然，稳定性条件(15)与Aboudi[1]的结果是一致的。

横向同性

如果是≤p，否则。

岩石圈结构和深部作用

其中，

岩石圈结构和深部作用

其中A，N，L，F和c是弹性常数。

同样，我们可以得到非均匀和其他特殊各向异性介质(如立方各向异性、正交各向异性等)的稳定性条件。

4三维各向异性情况下差分格式的稳定性条件

下面的稳定性条件可以通过前面二维情况下的分析获得:

定理2如果

Max [K1 ‖ A ‖+K2 ‖ E ‖+K2 ‖ Q ‖，f2(θ2，θ2，θ3)]≤1，则差分格式(3其中，函数f2(θ1，θ2，θ3)定义为:

岩石圈结构和深部作用

其中:、、、、、；

θ1，θ2和θ3满足:

岩石圈结构和深部作用

其中:a = 4k 1‖a‖b = 4k 2‖e‖d = 4k 3‖q‖g = 4k 4‖b+d‖e =

对于三维非均质介质的情况，可以用“冻结系数”法类似地分析稳定性条件。

5总结和讨论

有限差分法是数值模拟地震波传播的重要工具，其稳定性条件是提高计算速度的关键之一。然而，对于一般的二维和三维各向异性介质，对它们的差分格式和稳定性条件却少有系统深入的研究。本文给出了一种占用内存较少的快速差分格式，并系统地分析和推导了差分格式在一般均匀和非均匀各向异性下的稳定性条件。我们认为，本文的结果有助于各向异性数值模拟的发展，并为有限差分法的广泛应用提供了理论基础。

参考

[1]J .阿布迪.震源的数值模拟。地球物理学，1971，36，810~821。

[2]A.Bayliss，K.E.Jordan，B.J.LeMesurier和E.Turkel .计算弹性波的四阶精确有限差分格式. Bull.Seis.Soc.Am，1986，76，1115 ~ 1132 .

[3]D.C.Booth和S.Crampin .各向异性反射率技术:理论.地球物理学杂志，astr SOC，1983(a)，72，755~756 .

[4]D.C.Booth和S.Crampin .各向异性反射率技术:来自各向异性上地幔的异常到达波.地球物理学报，美国地球遥感学会，1983(b)，72，767~782 .

[5]V.Cerveny .非均匀各向异性介质中的地震射线和射线强度.地球物理学报. astr SOC，1972，29.1~13 .

[6]V.Cerveny和P.Firbas .非均匀各向异性介质中地震体波走时的数值模拟和反演.地球物理学报.美国地球物理学会，1984，76，41~51 .

[7]陈光海.各向异性非均匀介质中弹性波传播的数值模型:有限元方法.第54届SEG年会扩展摘要，美国休斯顿，1984，631~632。

[8]E.L.Faria和P.L.Stoffa .横向各向同性介质中的有限差分模拟。地球物理学，1994，59，282~289。

[9]H.Igel，P.Mora和B.Riollet .通过有限差分网格的各向异性波传播。地球物理，1995，60，1203~1216。

[10]K.R.Kelly，R.W.Ward，S.Treitel和R.M.Alford .合成地震图:有限差分法。地球物理学，1976，41，2~27。

[11]B.L.N.Kennett .弹性介质的理论反射地震图.地球物理学报，1979，27，301~321 .

[12]B.L.N.Kennett .分层介质中的地震波传播。剑桥:剑桥大学出版社，1983。

[13]D.Kosloff和E.Baysal .通过傅立叶方法的正演模拟。地球物理，1989，47，1402~1412。

陆，顾丽珍，陈。偏微分方程的差分方法。北京:高等教育出版社，1988。

[15]R.D.Richtmyer和K.W.Morton .初值问题的差分方法。纽约:Interscience，1967。

[16]C.Tsingas，A.Vafidis和E.R.Kanasewich .用有限差分法计算横观各向同性介质中的弹性波传播.地球物理学报，1990，38，933~949。

[17]J.Virieus.P-SV波在非均匀介质中的传播:速度-应力有限差分法。地球物理学，1986，51，889~901。

[18]张志军，何庆德，滕军伟，等..用有限差分法模拟二维横观各向同性介质中的三分量地震记录.勘探地球物理学报，1993，29，51~56 .

附录:引理

引理1是实系数方程λ 2-2dλ+1 = 0的充要条件，| d | ≤ 1是其根λ满足∣ λ| ≤ 1。

关于引理1的证明，见参考文献[14]。

引理2如果A∈Rn×n且A = A’，I是单位矩阵，则对于下面的方程。

∣λ2I—2Aλ+I∣=0，

‖A‖≤1是方程的根λ满足∣ λ | ≤ 1的充要条件。

证据:充足

由于a是实对称矩阵，T-1和T∈Rn×n的存在使得

t-1AT =诊断(d1，d2，…，dn)

根据引理2中的等式和上面的等式，我们可以得到

(λ2—2d 1λ+1)…(λ2—2dnλ+1)= 0

‖A‖≤1且ρ(a)≤a‖，▎▎ (a) ≤ 1。

根据，有

因此，对于满足方程的任何λ，都有∣λ∣≤1。

必要条件:由于∣λ∣≤1和a是实对称矩阵，我们可以从引理1得到:

∣di∣≤1,i=1,2,…,n

即ρ(A)≤1。

又因为A是正矩阵，ρ (a) = ‖ A ‖，即‖A‖≤1。

引理3如果A∈Rn×n且A = A’，则

(I)如果‖ i-a ‖≤ 1，则‖A‖≤2；

(ii)如果‖A‖≥0，则‖A‖≤2是‖ i-a ‖≤ 1的充要条件。

证明了:(I)因为A是实对称矩阵，所以存在T-1和T∈Rn×n使得T-1AT=diag(d1，d2，…，DN) ≡ D。

根据‖ I-A ‖≤ 1，有‖ T (I-D) T-1 ‖≤ 1。

很明显，I-D是正规矩阵，所以‖t(I-D)t-1‖=‖I-D‖≤1。

也就是

所以，也就是‖D‖≤2。

因为a是正规矩阵，所以

‖A‖= TDT-1‖= D‖，即‖A‖≤2。

(ii)在(I)中已经证明了必要条件，下面证明了充分条件。

根据(I)的证明过程:

‖A‖= D‖

因为‖A‖≤2且‖A‖≥0，所以有0 ‖ D ‖≤ 2，即最大∣ Di | ≤ 2。

这样我们就可以。

即:I-D‖= T-1(I-A)T‖≤1。

从A的对称性可以得到‖ I-A ‖≤ 1。