有限差分方法的论文文献

本文用与频率有关的时域有限差分法((FD)~2TD法),分析了在考虑人体色散特性时,脉冲电磁场与人体的作用。根据人体组织复介电常数实验原始数据在10MHz～3GHz的特性,提出一种2类组织(脂肪/骨/软骨、肌肉)简化色散模型,并与用12类组织色散模型的结果比较,验证了在分析人体内感生电流和比吸收量(SA)时,简化模型同样有效。计算结果显示,持续时间小于85ns、电场强度最大值达65kV/m的核爆炸典型脉冲作用下,人体内最大比吸收量(15ml/kg)小于ANSI/IEEE 标准限值()约三个量级。

杨顶辉

（石油大学地球科学系，北京100083）

滕吉文张中杰

（中国科学院地球物理所，北京100101）

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

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

1引言

地震波传播的数值模拟在地球科学中具有重要的意义。在各向异性地震模拟的各种方法中，基于Kennett研究工作[11,12]的反射率方法是最流行的数值技术之一。基于走时方程渐近解的射线追踪方法是模拟地震各向异性的另一种有效方法[5,6]。Kosloff等[13]利用Fourier方法模拟了地震各向异性。而Chen[7]则使用有限元方法模拟非均匀各向异性问题。虽然有限差分方法已被广泛应用于各向同性介质中的弹性波模拟，但是利用这种方法来模拟地震各向异性问题并不普遍。Tsingas等[16]利用有限差分算子发展了一种模拟算法以求解横向各向同性介质中的偏微分方程，这种算法是基于MacCormack型的分离格式[2]。Faria等[8]基于交错网格格式[17]，利用有限差分算法模拟了二维横向各向同性问题。最近，Igel等[9]基于褶积算法给出了一种模拟一般地震各向异性的有限差分算法。

快速和少存贮量是有限差分方法的优点。随着大尺度波场模拟的需要和大规模并行计算的发展，复杂介质或高维（二维和三维）模型的有限差分地震模拟不可避免。虚谱法因其空间算子准确地达到Nyquist频率而深受欢迎，然而虚谱法需要富利叶变换，从而对于三维各向异性模拟非常耗时。同时，采用富利叶变换意味着每一个网格点与其它的点相互影响。在某种意义上，这与动力学的局部弹性性质不一致。因此当我们为地震模拟设计有限差分格式时，考虑差分算子的局部性是必要的。另一方面，考虑差分算子的局部性对于提高算法的并行性非常重要。因为最邻近网格点间的信息交换最快，因而对于各向异性大尺度模型波场的模拟是可行的。

基于上述原因，本文针对一般各向异性介质中的地震波传播问题，给出了一种快速且占用内存少的有限差分方法。事实上，这是算法[10,18]的一种推广。

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

2有限差分格式

三维各向异性

各向异性介质中的运动平衡方程可成如下形式：

岩石圈构造和深部作用

其中：ρ是介质密度；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表示时间步长，并且，

岩石圈构造和深部作用

二维各向异性

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

岩石圈构造和深部作用

并且有下列差分格式：

岩石圈构造和深部作用

3稳定性条件

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

均匀介质情况下格式（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。

非均匀情况下的稳定性条件

在非均匀情况下，对于差分格式（5）的稳定性条件是难以确定的。然而根据偏微分方程的数值方法理论[15]，我们可以采用“冻结系数”法[14]来分析其稳定性。进一步地，我们给出非均匀介质情况下差分格式（5）的稳定性条件。

事实上，如果介质参数函数为连续有界的，则通常我们可以近似地将一个小的计算区域看成是均匀的，那么差分格式（5）可以退化成格式（6），这样在小区域范围内，我们可以像均匀情况一样获得稳定性条件。进一步根据介质参数函数的连续有界特性，我们可以获得差分格式（5）的稳定性条件为：

如果 那么

岩石圈构造和深部作用

否则，

岩石圈构造和深部作用

其中H（α，β）、α和β为：

岩石圈构造和深部作用

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

显然，通过计算格式（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=4k1·‖A‖，b=4k2·‖E‖，d=4k3·‖Q‖，g=4k4·‖B+D‖，e＝4k5·‖C+G‖，f=4k6·‖F+H‖。

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

5总结和讨论

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

参考文献

[1] simulation of seismic ～821.

[2] and fourth-order accurate finite-difference scheme for the computation of elastic ～1132.

[3] and anisotropic reflectivity technique:（a）,72,755～756.

[4] and anisotropic reflectivity technique:anomalous arrivels from an anisotropic upper （b）,72,767～782.

[5] rays and ray intensities in inhomogeneous anisotropic ～13.

[6] and modelling and inversion of travel-times of seismic body waves in inhomogeneous anisotropic ～51.

[7] numerical model of elastic wave in anisotropic inhomogeneous media:finite element 54th SEG Annual Meeting Expanded Abstracts,Houston,～632.

[8] and modelling in transversely isotropic ～289.

[9] and wave propagation through finite-difference ～1216.

[10] and seismograms:A finite-difference ～27.

[11] reflection seismograms for an elastic ～321.

[12] wave propagation in stratified University Press,1983.

[13] and modelling by a Fourier ～1412.

[14]陆金甫，顾丽珍，陈景良.偏微分方程差分方法.北京：高等教育出版社，1988.

[15] and methods for initial value York:Interscience,1967.

[16] and wave propagation in transversely isotropic media using finite ～949.

[17] wave propagation in heterogeneous media:Velocity-stress finite-difference ～901.

[18] al..Simulation of 3-component seismic records in a 2-dimensional transversely isotropic media with ～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=diag（d1,d2，…，dn）

根据引理2中的方程及上述等式，可获得

（λ2—2d1λ＋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，即max∣di｜≤2

所以我们可获得 。

即：‖I—D‖＝‖T-1（I—A）T‖≤1

由A的对称性可得‖I—A‖≤1。