复杂地表地质条件下地震波数值模拟综述.pdf

复杂地表地质条件下地震波数值模拟综述 王雪秋1, 2,孙建国1, 2,张文志3 1. 吉林大学 地球探测科学与技术学院,吉林 长春 130026; 2.国土资源部 应用地球物理综合解释理论开放实验室 波动理论与成像技术实验室,吉林 长春 130026; 3.中国国土资源航空物探遥感中心,北京 100083 摘要通过对复杂地质地表条件的分析研究,总结了近年来在这方面地震波场数值模拟的新方法。复 杂地表的类型很多,对地震资料的影响程度不同,要想完全描述近地表的影响因素是不可能的。文中所介 绍的模拟方法在处理这些问题上都是各有侧重、 各有优势,目的是寻找合适的方法解决适合的问题。 关键词数值模拟;映射;不规则网格;起伏地表 中图分类号 P315. 31 文献标识码 A 作者简介王雪秋19742 , 女,辽宁沈阳人,讲师,博士研究生,主要从事地震波理论及波场正演模拟等方面的研究. The State-of the-Art in NumericalM odeling Including Surface Topography WAN G Xue2qiu 1, 2, SUN Jian2guo1, 2, ZHAN GW EN 2zhi3 1. College of GeoExp loration S cience and T echnology,J ilin U niversity,Changchun130026; 2.L aboratory of W ave T heory and Im ag ing T echnology,Open R esearch L aboratory f or Integ rated Geophysical Interp retation T heory of theM inistry of L and and R esources,Changchun 130026,China; 3.China A ero Geophysical S urvey and R em ote S ensing Center f or L and and R esourse,B eijing100083,China Abstract A s a tool for understanding the regulation if w ave propagation in complex media and exam ining the effects of different kinds of s, seism ic numericalmodeling plays an i mportant role in modern geophysical prospecting.This paper briefly introduced several modeling s including surface topography in recent years .Our ai mis to find the suitable to solve the adapted questions . Key words numericalmodeling; mapping; irregular grids; surface topography 0 引言 消除地表地形对地震资料的影响一直都是地震 勘探中的一个世界性的难以克服的问题,尤其是在 中国。众所周知,中国的地形地势是非常复杂的,但 是在当今对油气资源大力需求的形式下,地震勘探 的前缘正在逐步向复杂地表地区转移,尤其是在西 部的沙漠、 高原、 山地和丘陵地区,地震工作已经大 面积的展开了。所以对这些地区的地表地形影响的 研究就显得日益迫切。 一般情况下,对地表地形影响 的消除重要是用静校正方法。但是这里有两个问题 需要解决其一是静校正量的确定,对于复杂的地表 条件,这是一个很难确定的问题,需要做大量的正演 模拟和野外调查;另外一个因素是致命的,目前的静 校正方法是以几何地震学的假设前提地表水平,表 层介质均匀为前提的 , 显然复杂地表地形的影响是 不在这个假设之内的[1 ～3]。 近年来,地震勘探所涉及的复杂地表包括近地 表的不均匀性引起近地表散射 , 起伏地表、 厚的黄 土覆盖层,地形起伏变化复杂,地表高差大山地。 山地的特点是地形剧烈起伏,岩石出露,山间地表松 软,低速层厚度变化大,山前戈壁砾石区近地表速度 变化大,高速岩层直接出露地表不存在低速带,速 度横向变化大,近地表射线传播路径不垂直复杂地 表使速度的提取精度受到影响,另外复杂的地表地 形散射也影响了深层信息的提取。例如世界上公认 第35卷 专辑吉 林 大 学 学 报地 球 科 学 版Vol135 Sup1 2005年7月Journal of Jilin U niversityEarth Science EditionJuly 2005 的地震勘探困难地区黄土高原,主要表现在两个方 面巨厚的黄土,干燥疏松,地震波吸收衰减严重,激 发和接收条件差;表层结构复杂,地形高差大,黄土 层厚度和低速带厚度、 速度横向变化剧烈,静校正问 题突出,测线布设以沿沟弯线为主。黄土塬地区,一 方面是激发能量弱,另一方面是黄土对地震波的吸 收衰减严重。 相对于白垩系砂泥岩而言,黄土层对地 震波的吸收是严重的。 比较同一检波点,接收在基岩 与黄土中激发折射波初至的能量,含膏砂泥岩是黄 土的5～15倍,相差14～24 dB,地震波在黄土层中 的衰减是基岩的2～8倍。结论是在黄土上激发区无 论什么埋置条件,资料总是很差,红土的激发能量比 黄土强。 随着油气勘探难度的加大如复杂构造、 复杂岩 性、 复杂地表条件 , 传统的基于平缓构造勘探的地 震数据采集、 处理和解释方法技术面临一系列挑战, 许多技术在这类复杂地区已经不再适用。地震波传 播数值模拟技术是研究地震波传播规律的有效途 径,通过地震波传播数值模拟分析,弄清地震波在复 杂地质构造中的传播特征,可以根据不同地区的实 际情况,采用更加合理的野外观测系统如测线方向 和排列长度 , 使用更加合理的处理流程和处理参 数,进行合理的地震解释。 1 复杂地表模拟定性分析 复杂地表的类型很多,对地震资料的影响程度 也不尽相同,要想完全描述近地表的影响因素是不 可能的。为了模拟地表严重起伏和表层强非均匀性 对地震波传播影响,算法的选择要符合以下几点 1不平的自由地表面必须在几何上精确描述并且 自由地表面容易用算法实现 ; 2 所有前向、 反向散 射,包括在局部复杂速度构造层内和强不均匀介质 内的鸣震,必须在模型正演过程中得到正确处理。 对于地球中地震波的传播来说,地震数值模拟 是一项重要的技术。其目标就是在假定地质构造已 知的情况下,预测能够在地面接收到的理论地震图。 对于地震解释和地震反演方法来说,这是一种非常 重要的技术。其另一个比较重要的方面是为地震勘 探提供评估和设计。模拟方法有很多种,可以分成3 类直接法、 积分方程法和射线追踪法。 直接法就是解波动方程,即将地质模型近似放 在一个数值网格上,也就是一个有限数字离散点的 模型。 这些方法也可以叫做网格法或全波动方程法, 后者的叫法是因为方程的解给出了全波场的特征。 直接法对于介质的的变化没有什么限制,并且在足 够小的网格上能够得到非常精确的解。这种方法能 够处理含不同流体的问题,而且非常容易得到波场 传播快照辐射图 , 这对结果的解释非常重要。 其缺 陷就是相对于其他方法来说,耗机时较多。 积分方程法是基于点源波场波级数的积分表 达。 该方法起源于惠更斯原理,详细研究惠更斯原理 可以认为波场一方面可以表示为一系列点源所形 成的波场的叠加,另一方面可以表示成由边界点源 所形成的波场的叠加。 惠更斯原理现在仍然在使用, 称这两种方式为体积分方程和边界积分方程,并且 各有自己的适用范围。这两种方法在应用上比之于 直接方法有更多的限制。 然而对于特定的几何形态, 例如均匀介质中包含异常绕射体问题、 井间问题及 微裂隙等问题,积分方程法的计算效率和精度都很 高。 由于具有良好的解析性,该方法对于建立在 Born近似基础之上的成像方法也很有用。 渐进法或射线追踪法是地震模拟和成像领域中 应用最广泛的技术。 由于没有考虑全波场特征,所以 是一种近似方法。 但是从另外一个角度看,这种方法 也是本文所讨论的方法中效率最高的方法。 综上所述的方法中,能够适应复杂地表地质条 件的模拟方法一般选取直接法波动方程方法。这 是由实际问题的要求和该方法本身的特点所决定 的[4 ～6]。 2 正演模拟定量研究 波动方程的解法有很多,总体上以网格法为主, 包括有限差分、 有限元、 傅立叶变换比等方法及这三 类方法的联合方法等。其基本特点是能够提供的波 场快照图可以帮助研究者细致分析波场在各种复杂 介质内部的反射、 透射、 绕射、 散射以及能量的衰减 等运动学和动力学的各种细节特征。 由于实际工作中所模拟的介质不同,所用的模 拟方程也不一样。 根据模拟方程的不同,波动方程数 值模拟主要有声波模拟、 弹性波模拟、 粘弹性波模 拟以及裂隙和孔隙弹性的iot模拟等。 波动方程模拟 方法的数学实现是数值模拟的一项重要研究内容, 它需要极高的计算机资源,这也是长期制约该方法 走向生产实际的最主要因素,并因此产生了针对各 种不同地质问题的数值实现方法。例如有限元 FE、 有限差分FD、Fourier变换、 边界元BE、 31专 辑 王雪秋,孙建国,张文志复杂地表地质条件下地震波数值模拟综述 谱元SEM、 广义屏等方法。这些方法的优劣见表 1。 表1 各种波动方程解法比较 Table 1 The comparison of s for wave equation 名称优点缺陷 有限元法处理地表起伏问题强计算量大,计算效率低 傅氏变换法运算速度快难应付横向变速 谱元法处理边界灵活计算量大,计算效率低 边界元法表征问题准确,适应各种 边界问题 构造方程系非常大,计算 量大 广义屏法运算速度快介质速度变化小和入射 角小 普通有限差 分法 运算速度快,实现简单频散严重、 难应付起伏地 表问题 高阶有限差 分法 精度高,运算速度较快, 实现简单 难应付起伏地表问题 这些解法的基础都是在离散网格上计算基于各 种方法实现的波动方程代数方程或方程组。因此解 决复杂地表问题的上,基于网格化这个层次上可以 分为两类一类是从模型出发,改变方程,建立适应 这种情况的方程,如映射网格,可调节网格,坐标变 换等;另外一种方法就是建立适当的边界条件。 在处理起伏界面问题时,离散化数值计算只能 将连续界面变成阶梯状界面,在大间距的网格上做 各种数值计算必然会产生不必要的人为绕射和由于 界面条件不准引起的积累误差。为了消除这些绕射 和误差,产生了许多方法[7]。 211 映射网格 将直网格改成曲网格是一种很自然的想法。 在曲网格上做数值计算的思想最早由Fornberg 图1 直网格变曲网格 Fig11 Mapping of a rectangular grid onto a curved 提出,但是Fornberg的映射技术只能使网格沿着主 要的界面变化。N ielsen等采用了计算流体力学中常 用的分块映射技术使网格可以沿着所有界面变化, 提高了算法的精确性。使数值网格沿着所有地下界 面变化是N ielsen的主要思想。为实现这一目标,可 将实际的几何构造影射成某一曲坐标系下的规则图 形,比如矩形。再利用插值法产生区域内部的各个 点,从而生成所需的曲网格图 1 [1, 7]。具体实现如 下将直角坐标系z-x下的任意一个不规则多边形 的边分成4段,每段对应于四边形中的一条边。将2 条相对的边分别映射成曲坐标系N-G域下的2条平 行于G轴的对边,类似地将另2条相对边映射成曲坐 标系下的2条平行于N轴的对边。通过这4条边界, 任意一个直角坐标系的不规则区域可映射成了一个 矩形,相应地建立了一个曲坐标系,然后在曲坐标系 下建立一个规则网格,就可以在规则网格上进行数 值计算了。一般来说网线常常取成整数,如N 1, 2, ⋯,G 1, 2,⋯,曲坐标系下网点的坐标为N,G,而其 实际坐标值为xxN,G , zzN,G。 图形内部点实 际坐标值采用一种多维插值法超限插值法来计 算,即利用4条边界上的点产生内部点的一种插值 方法。 r _ N,G ∑ 2 n 1 UnNr _ Nn,G ∑ 2 m 1 WmGr _ N,Gm ∑ 2 n 1∑ 2 m 1 UnNWmGr _ Nn,Gm , 式中r _ x e _ xze _ z,e _ x和e _ z是直角坐标系下的 单位向量;Un和Wm是插值函数。为方便起见,我们取 线性插值函数 U1N N2-N N2-N1, U2N N-N1 N2-N1, W1G G2-G G2-G1, W2G G-G1 G2-G1。 为了使网格沿所有的地下界面变化,并防止因 41 吉 林 大 学 学 报地 球 科 学 版 第35卷 界面过于复杂而造成的网格过于不规则或过于倾 斜,我们采用了在计算流体力学中应用广泛的分块 方法。 每一层被分成若干相邻的块,将每一块单独映 射生成网格后再把所有的块连起来形成统一的整 体。连接时应注意网格的连续性。 因为数值运算是在曲坐标系下的规则网格上进 行的,所以应把直角坐标系的声波方程变换成曲坐 标系下的方程。直角坐标系下的声波方程为 Q5 u 5t - 5p 5x sxx,z , Q5 w 5t - 5p 5z szx,z , 5p 5t -K5 u 5x -K5 w 5z 。 式中p -R表示压力R是应力 ; u和w表示位移 向量的分量;Q表示密度;KQv2v表示波的传播速 度表示L ame常数;sx和sz表示震源。将x-z坐标 系下的声波方程转换到曲坐标系N-G下有 Q5 u 5t -Nx 5p 5N- Gx 5p 5G sxN,G , Q5 w 5t -Nz 5p 5N- Gz 5p 5G szN,G , 5p 5t -K Nx 5u 5N Gx 5u 5G Nx 5w 5G Gz 5w 5G , sxN,G Nx 5f 5N Gx 5f 5G ht, szN,G Nz 5f 5N Gz 5f 5G ht。 系数可由隐函数求导法则导出 Nx zG J ,Nz xG J ,Gx zN J ,Gz xN J 。 这里J是Jacobian行列式 JxNzG-xGzN。 对插值函数r _ N,G分别关于N和G求导数,则可 得xN,xG,zN和zG。 5 5Nr _ N,G ∑ 2 n 1 5 5N UnNr _ Nn,G ∑ 2 m 1 5 5N WmGr _ N,Gm - ∑ 2 n 1∑ 2 m 1 5 5N UnNWmGr _ Nn,Gm , 5 5Gr _ N,G ∑ 2 n 1 5 5G UnNr _ Nn,G ∑ 2 m 1 5 5G WmGr _ N,Gm - ∑ 2 n 1∑ 2 m 1 5 5G UnNWmGr _ Nn,Gm。 这样,映射后在曲坐标系下得网格是一个规则 网格,波动方程就可以在这样的网格上进行计算。 212 不规则网格 非规则网格差分算子与上面的分段映射方法有 相似之处,通过在离散混合变量弹性波方程时引入 非规则网格差分算子,应用于应力-速度方程可以 适用于任意形状的边界和内部边界[3, 8, 9]。该方法的 本质就是一个精确处任意形状自由表面的边界条件 的坐标变换问题。 这个方法从任意四边形网格着手, 在非正交网格上定义任意四边形网格差分算子,由 任意四边形的4个角点的函数值计算函数在域内的 空间导数,若将x,y坐标系下的任意四边形映射 到N,G坐标系下的标准四边形,则x,y坐标系 下任意四边形网格内的差分微分计算可以转化为 在N,G坐系下的正方形网格上进行,定义映射函数 为 x 1 4 1 N 1 Gx1 1 4 1 - N 1 Gx2 1 4 1 - N 1 - Gx3 1 4 1 N 1 - Gx4, x 1 4 1 N 1 Gz1 1 4 1 - N 1 Gz2 1 4 1 - N 1 - Gz3 1 4 1 N 1 - Gz4。 其中xi,zi i 1, 2, 3, 4是x,y坐标系下四边形4 个角点的坐标。该映射函数可以将图中左侧的任意 四边形映射为右侧的边长为2的标准正方形,若在 N,G坐标系下对函数Px,z求导,则有 5P 5N 5P 5x 5x 5N 5P 5z 5z 5N , 5P 5 5P 5x 5x 5N 5P 5z 5z 5N 。 由上式可以得到 5P 5x 5P 5z 5x 5N 5z 5N 5x 5G 5z 5G - 1 5P 5N 5P 5G 。 设函数Px,z在四边形4个角点处的函数值Pi i 1, 2, 3, 4,则N,G坐标系下函数Px,z可以 51专 辑 王雪秋,孙建国,张文志复杂地表地质条件下地震波数值模拟综述 用如下双线性插值描述 P 1 4 1 N 1 GP1 1 4 1 - N 1 GP2 1 4 1 - N 1 -GP3 1 4 1 N 1 - GP4。 对上式求导,得 5P 5N 5P 5G 1 4 1 G-1 -G-1 G1 -G 1 N1 -N-1 N-1 -N P1 P2 P3 P4 。 将这个式子代入上面的矩阵中,即可得到任意四边 形4个角点的函数值,计算函数在域内空间上导数 的差分算式。 类似的方法还有尧德中等利用六角形网格的周 期采样,通过线性变换直接利用傅立叶变换方法求 解波动方程,也可以得到处理复杂边界条件的目 的[10]。 另外由M artin Kaser等人提出的不规则三角网 格,可以任意剖分不规则边界,通过一个不规则加权 函数来离散波动方程,来求解复杂条件下的波动方 程的解,包括复杂地表条件[11]。 其实改变网格形状或大小的不规则网格还有很 多种,例如最小网格、 不均匀网格、 加权网格等很多 Dong2Joo M in,Changsoo Shin,Byung2Doo Kwon and Seunghw an Chung,2000; Ivo Oprsal andJiliZahradnik,1999;Sophie2A delaide M agnier,peter M ora, and A lbertTarantola, 1994; Jiri A ahradnic, Patrick O ‘ leary, and James Sochackl,1994; Zhang Jianfeng and L iu T ielin, 1999; Zhang Jianfeng , 1999; 。在有些情况下,这样 的剖分方式更有助于算法的实现和精度的保 证[12 ～15]。 213 可调节网格 这种方法的思路与正常的差分矩形网格一样, 但是网格的横向或纵向间距是可以调节的,也就是 说在不同的界面上可以调节网格的大小来适应地形 或构造的变化,通过一个可以调节的加权因子来控 制网格的大小 5 5z A5 f 5x 5 5z g, 5 5z g 0, 0≈ 2 g0, 1 2 -g0, 1 2 D Z - 1 D Z0。 这实际上是通过对节点上的参数进行加密来实现各 种复杂构造及地表条件的模拟。 214 坐标变换 这种方法实际上是对相对比较简单的线性斜坡 式的地表进行的一种比较简单的变化,可以通过图2 来实现。 图2 斜坡坐标变换 Fig12 A ramp with slope to illustrate algorithm 215 不同阶差分法 用波动方程有限差分方法做正演模拟时,如果 使用二阶差分格式,网格间距必须取得很小,以保证 计算精度和稳定性,而这样做明显地降低了计算效 率。 为此,Dalain和M ufti提出了高阶有限差分方 法,可以适用于复杂介质强横向变速情况下的正演 模拟。在高阶方程的情况下,网格间距可以取得很 大,但计算精度并不比二阶差分方程小网格间距时 低。而网格间距增大可以显著提高计算效率。一般 称截断误差高于四阶的差分方程为高阶差分方程。 三维声波方程 5 2u 5x 2 5 2u 5y 2 5 2u 5z 2 1 v 2x ,y,z 5 2u 5t2 , 其中ux,y,z,t为地表记录的压力波场,vx,y,z 为纵横向可变化的介质速度。在局部速度变化很缓 或不变时,截断误差为Ox M ,y M ,Z M ,t4的统 一的三维正演模拟高阶差分方程是 Ltt 2ut - ut-t vt 25 2u 5x 2 5 2u 5y 2 5 2u 5z 2 2 4 vt 45 4u 5x 4 5 4u 5y 4 5 4u 5z 4 4 4 rt 4 5 4u 5x 25y2 5 4u 5y 25z2 5 4u 5z 25x2 Ot4。 这是一个推导用于正演模拟高阶差分方程的起始方 程。 它在时间方向上截断误差为Ot4 , 而空间的截 61 吉 林 大 学 学 报地 球 科 学 版 第35卷 断误差则要根据需要来定至少是四阶以上的。另 外,可以根据需要来组合不同阶次的差分格式。 当局 部速度横向变化很大,即vvx,y,z时,截断误差 为Ox M ,y M ,Z M ,t4的统一的三维正演模拟高 阶差分方程为 Ltt 2ut - ut-t vt 2 5 2u 5x 2 5 2u 5y 2 5 2u 5z 2 t2 6 5v 5x 2 v 5 2v 5x 2 5v 5y 2 v 5 2v 5y 2 5v 5z 2 v 5 2u 5z 2 5 2u 5x 2 5 2u 5y 2 5 2u 5z 2 v 3 t4 3 5v 5x 5 3u 5x 3 5 3u 5y 25x 5 3u 5z 25x 5v 5y 5 3u 5x 25y 5 3u 5y 3 5 3u 5z 25y 5v 5z 5 3u 5x 25z 5 3u 5y 25z 5 3u 5z 3 1 6 vt 45 4u 5x 25y2 5 4u 5y 25z2 5 4u 5z 25x2 1 12 vt 45 4u 5x 4 5 4u 5y 4 5 4u 5z 4 Ot4。 其中十阶差分方程的系数是当M 10时, X0 - 518544445,X1 31333333, X2 - 014761901,X3 0107936513, X4 - 01009920621,X5 010006349185。 216 边界条件 总结起来网格法的边界条件有四种高衰减区 法、 傍轴近似法、 吸收边界法、 多次入射法等,但是在 起伏地表时,上界面不再是平界面,上述的方法在理 论上和实现起来都有一定的困难。但是陈伟提出的 将边界条件以波动方程的形式给出,在一个地表起 伏的计算区域内,边界方程写成带特征值的特征矩 阵形式 5W 5t L - 1 AM LA 5W 5x L - 1 CPLC 5W 5z 。 这里L是特征矩阵。 根据边界起伏特点,处理这 个方程可以得到不同间断处的边界条件,这样有限 差分解波动方程时的边界问题就转化为在边界处解 同样形式的波动方程的有限差分问题。 217 频散处理 在地形起伏情况下,为减少地形和表层不规则 性对体波和面波的散射,需要更有效的减少数值频 散,提高模拟精度。 数值频散是由于数值离散而产生 的,它使得具有不同频率的地震波表现为不同的相 速度,从而使地震波发生弥散,影响数值模拟的效 果。 数值模拟和差分精度存在直接的关系,随着差分 精度的提高,数值频散逐渐减小。在波动方程计算 中,尽管数值频散不可避免,但是通过采用高阶差分 解法可以提高计算精度,并可以尽量减少数值频散。 在声波方程中,波场的空间微分的高阶差分近 似为 5 2f 5x 2 1 x 2∑ M m 1C m m[fxmt - 2fx fx-mt ] 。 式中,高阶差分系数C m m可以根据不同的差分精度 求得。 研究发现,地震波相速度c和波数k之间具有如 下关系 c2 c20 ≈1 2- 1 M U 2M 2 M 2∑ M m 1C m mm 2M 2 通过上式的分析可以发现,只要满足Ukx 2P K, x≤1,即x≤ K 2K实际上, x还可以大一些 , 也就 是说只要空间离散时一个波长包含几个x,随着差 分精度 2 M 的提高,上述高阶差分解法产生的数值 频散就会逐渐减小。 3 结论 对于复杂的地表条件地震勘探来说,消除地表 的影响一直是个难以克服的棘手问题。所以解决问 题的第一步应该是通过数值模拟来定量分析,然后 才能根据实际情况采取相应的接收处理手段。以上 介绍的方法都是这个领域寻近10年公开发表的研 究方法。我们的目的是寻找合适的方法解决适合的 问题。 到目前为止,关于复杂地表的地震波场数值模 拟技术尽管有了很大的发展,但是相对来说,还是处 在一个比较新的研究阶段,随着石油勘探技术的不 断发展,这项研究将是今后一个时期内的热点研究 课题。 参考文献 [ 1] 赵景霞,张叔伦,孙沛勇.曲网络伪谱法二维声波模 拟.石油物探, 2003, 421 125. [2] 裴正林.任意起伏地表弹性波方程交错网格高阶有限 差分法数值模拟.石油地球物理勘探, 2004, 39 6 6292634. [3] 张剑峰.弹性波数值模拟的非规则网格差分法.地球 物理学报, 1998, 41增刊 3572366. [ 4 ] Claudia Kerner.M odeling of soft.sedi ments and liquid2solid interface modified wavenumber summa2 71专 辑 王雪秋,孙建国,张文志复杂地表地质条件下地震波数值模拟综述 tionandapplication[ J ].Geophysical Prospecting, 1990, 38 1112137. [5] Kelly K R,W ard R W , Sven T reitel, et al .Synthetic seismogram s afinitedifference approach[ J ]. Geophysics, 1946, 411 2227. [6] L iu H sui2lin, M ichael Prange, Francois Daube.Fini2 tedifference modeling of sonic wavefields w ith a dipping interface [J ].Geophysics,1997,62 4 127021277. [7 ] N ielsen P, If F, Berg P, et al .U sing the pseudos2 pectral technique on curved grids for 2D acoustic forward modeling [ J ].Geophysical prospecting, 1994, 423 3212342. [8] Zhang Jianfeng, L iu T ielin.P2SV2wave propagation in heterogeneousmedia grid [J ].Geophys J Int, 1999, 136 4312438. [9] Zhang Jianfeng. Q udrangle2grid velocity2stress finite difference for poroelastic wave equation [J ]. Geophys J Int, 1999, 139 1712182. [10] 尧德中,阮颖铮.六角形网格波动方程数值模拟德傅 氏变换算法[J ].计算物理, 1994, 111 1012106. [11] M artin Kaser, Heiner Igel . N umerical si mulation of 2D wave propagation on unstructured grids using explicitdifferentialpoerators [ J ].Geophysical Prospec2ting, 2001, 49 6072619. [12] Dong2Joo M in, Chang2soo Shin, Byung2Doo Kwon, et al . I mprovedfrequency2domainelastic wave modelingusingweighted2averagingdifference operators [J ].Geophysics, 2000, 653 8842895. [13 ] Ivo Oprsal, Jili Zahradnik.Elastic finite2difference for irregular grids [J ].Geophysics, 1999, 641 2402250. [14] KindelanM , KamelA , Sguazzero P. ShortNote On the construction and efficiency of staggered nume2 rical differentiators forthe wave equation[ J ]. Geophysics, 1990, 1551 1072110. [15 ] M ichael Korn.Computation of wavefields in verti2 cally inhomogeneous media by a frequency domain finite -difference and application to wave propagation in earth models w ith random velocity and density perturbations [J ]. Geophy J R astr Soc, 1987, 88 3452377. 81 吉 林 大 学 学 报地 球 科 学 版 第35卷