一阶弹-方程交错网格高阶差分解法f=生波稳定性研究

《一阶弹-方程交错网格高阶差分解法f=生波稳定性研究》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、第43卷第6期地球物理学报Vd43．No62000年11月CHINESEJOURNALOFGEOPHYSICSN．2O00[文章编号]00015733{2000)060856—09[中囤分类号]1'631一阶弹-f=生波方程交错网格高阶差分解法稳定性研究平耕董良国马在田曹景忠(同痹夭荸吾薄洋丽面赢验室．上海200092)7、／[摘要]稳定性问题是数值求解波动方程的基本问题．文中对三堆横向各向同性介质中一骱弹性波方程交错网格高阶差分解法的稳定性进行了分析．给出了不同精度差分方程统一的稳定性条件，证明了

2、三维TI介质中一阶弹性波方程交错网格高骱差分解法的稳定性由弹性波在x、y、z三个方向上的Courant数共同决定．最后通过几种精度差分方程的稳定性条[关键词]横嚣向鸳壁，茎塾星丝件，煎睑墓羞盐这，交错网格．的坚’l塑!!二一1己I吉协触在弹性波正演模拟中，除了使用二阶方程外，常常采用一阶速度一应力弹性波方程，其主要优点是勿需对弹性常数进行空间微分在具体的有限差分解法上，除了规则网格外，一种较为先进的交错网格最早由Madariaga，R．[11提出，Virieux在模拟各向同性介质中的SH波和PSV波

3、时也使用了这种差分网格，其差分精度为O(at+△)，在不增加计算工作量和存储空间的前提下，和常规差分网格相比局部精度提高了4倍，且收敛速度也较快J．Lev'anderI5又将这种差分网格的精度提高到O(△+△z)．为了进一步提高差分精度，减小网格频散，本文作者将速度(应力)对时间的奇数阶高阶导数转化为应力(速度)对空间的导数，从而在不增加所需内存量的前提下，将交错网格和高阶差分法有机结台，成功地运用到求解二维TI介质一阶速度一应力弹性波方程中，并在实际模型的弹性波模拟中取得了较好的模拟效果．这种求解

4、一阶弹性波方程的交错网格高阶差分解法，与常规的有限差分法以及交错网格低阶差分法相比，网格频散显著减小，精度明显提高．在保证一定的差分精度的基础上，可以取较大时间步长和空间网格间距，提高了计算效率．由于差分计算过程中数值参数选择不台理，可能产生无物理意义的按指数增大的数值计算结果J，造成模拟结果网格频散严重，影响对问题的分析，严重时会造成溢出而使计算f收稿日捐]1999—11—09收到，2I]0o一【l5—24收刊修定稿．【基盒项目]教育部重点科技基盎项目和薄洋863—820青年基盎项目【作者简介]董

5、良国，男．1966年生1990年毕业于同济大学勘探地球物理专业，现为该校在职博士研究生．主要从事地震渡的传播以厦地震数据灿理等方面的教学和科研工作辫i6期董良国等：一阶弹性波方程交错网格高阶差分解法稳定性研究857无法进行因此，对一种数值解法，需要知道计算稳定的离散参数区域，即分析解法的稳定性2一阶弹性波方程的高阶差分形式2．1三绻一阶弹性波方程在不考虑体力的情况下，由速度一应力表示的三维一阶弹性波方程可为：Qv()，(1)这里，U=(口，口，，，，，，，)，为速度，r为应力，Q是含有介质弹性常数以

6、及空间微分算子D、D和D的9×9阶方阵2．2一阶弹性波方程交错网格高阶差分形式将【，(tq-等)和【，(t-)通过Tay】0r公式展开后相减，可得2M阶时间差分精度的差分近似u等)一u等)=z薹M)'其中，△￡为时间步长．当M=1时，(2)式即为传统的二阶精度时间差分近似．为了提高时间差分精度，同时避免直接计算(2)式中的(f)时所涉及过多时间层而需求过多的内存量，利用方程(1)可以完全准确地将速度对时间的任意奇数阶高阶导数转为应力对空间的导数，将应力对时间的任意奇数阶高阶导数转为速度对空间的导数．

7、这样，在计算一个时间层上速度(应力)时，只需要前一时间层的速度(应力)场，以及中间节点上的应力(速度)场，不需要过多的时间层，从而节省内存．由(1)式可知氅1一I：Qz一【一，(、⋯)，(⋯3)为了提高空间差分精度，采用交错弼格法得到的2N阶空间精度的近似式为型3x=耋c{_，[zq-等(2n一1)]一，[一警(2一1)]}+0((4)其中，△z为差分网格间距，c’为不同阶的交错网格差分系数(该系数的确定方法可参见文献[6])．将(3)式中Q精确的空间微分算子DD和D用(4)式的2N阶空间精度的近似

8、微分算子代替后，变成矩阵R，这样(3)式就变为如下2N阶空间差分精度的差分形式≈R【，(￡)，(5)将(5)式代入方程(2)，即可以得到2M阶时间、2N阶空间差分精度的三维一阶弹性波方程的交错弼格高阶差分形式u譬)一u(卜等)=z耋R2"-1U㈤+地球物理学报43卷O(／,tM+△+△+△2)，(6)3高阶差分方程稳定性分析方法方程(6)两边对时间求导后，分别对时问和空间进行Fburier变换，可得[Esin2(~1)+](，∥)：0，(7)这里E为单位矩阵，G(，k