波动方程数值模拟的三种方法及对比

第２４卷第２期

２００９年４月（页码：３９１～３９７）

地球物理

ＩＮ

学进展

ｖ０１．２４，Ｎｏ．２

Ａｐｒ．２００９

ＰＲＯＧＲＥＳＳＧＥＯＰＨＹＳＩＣＳ

皮红梅，蒋先艺，刘财，等．波动方程数值模拟的三种方法及对比．地球物理学进展，２００９，２４（２）：３９ｌ～３９７，ＤｏＩ：１０．３９６９／ｊ．

ｉｓｓｎ．１００４—２９０３．２００９．０２．００４．ＰｉＨ（ｉｎ

Ｍ，ＪｉａｎｇＸＹ，ＬｉｕＣ，Ｐｆ口Ｚ．Ｔｈｒｅｅｍｅｔｈｏｄｓｏｆｗａｖｅｅｑｕａｔｉｏｎｆｏｒｗａｒｄｍｏｄｅｌｉｎｇａｎｄｃｏｍｐａｒｉ５０ｎｓ．ＰｒＤｇ化５ｓ锄ＧＰｏ户＾∥．

Ｃｈｉｎｅｓｅ），２００９，２４（２）：３９１～３９７，ＤｏＩ：ｌＯ．３９６９／ｊ．ｉｓｓｎ．１００４—２９０３。２００９．０２．００４．

波动方程数值模拟的三种方法及对比

皮红梅１’２，

蒋先艺２，

刘

财¨，

王成祥２，

姜绍辉２

（１－吉林大学地球探测科学与技术学院。长春１３００２６｝

２．中国石油东方地球物理公司物探技术研究中心，涿州０７２７５０）

摘要波动方程数值模拟方法是研究地震渡场传播的一种重要手段，本文采用交错网格高阶有限差分方法分别对

双程声波方程和双程弹性波方程进行了波场数值模拟，并且根据定位原理采用傅立叶有限差分算子进行了单程波方程数值模拟，在分析定位原理的基础上，对其计算过程稍作修改，将延拓到地面的波场直接由每个检波点接收，无需横向叠加过程，得到了单程声波方程共炮记录．基于不同波动方程的数值模拟结果表明，双程波方程结果包含直达波、多次波等干扰渡，信噪比低；单程波数值模拟结果只包含了介质分界面的一次反射波，信噪比高，但对于大角度入射渡误差较大，并且对于同一个地质模型而言，双程弹性波方程计算速度最慢，双程声波方程次之，单程声波方程计算速度最快．因此对于复杂地质模型，三种模拟方法可以取长补短，综合应用．关键词波动方程数值模拟，交错网格，高阶有限差分法，定位原理，单程波

ＤＯＩ：ｌＯ．３９６９／ｊ．ｉｓｓｎ．１００４—２９０３．２００９．０２．００４

中图分类号

Ｐ３１５

文献标识码

Ａ

ＴｈｒｅｅｍｅｔｈｏｄｓｏｆｗａＶｅｅｑｕａｔｉｏｎｆｏｒｗａｒｄｍｏｄｅｌｉｎｇａｎｄｃｏｍｐａｒｉｓｏｎｓ

ＰＩＨｏｎｇ—ｍｅｉｌ”，

（１．ＣＤｚＺＦｇＰ

ＪＩＡＮＧＸｉａｎ—ｙｉ２，

ｏ厂ＧＰ０８ｚｐＺｏｒ乜￡ｉｏ靠Ｓｃ｛Ｐ，ｌｃｅ

ＬＩＵＣａｉ¨，

ｎ行ｄ

ＷＡＮＧ

Ｃｈｅｎｇ—ｘｉａｎ９２，

ＪＩＡＮＧＳｈａｏ＿ｈｕｉ２

Ｔ■ｆ＿Ｉｌ雅ｏｚｏｇｙ，．，订ｉ以ＵｈｉｔＪＰ，百ｆ￡，，（冼口”ｇｃ｜Ｉｌ配种１３００２６，ｃ＾ｆ，ｌ口；

２．Ｒ＃ｓ∞”九＆ＤＰ御ＺＤ户ｍ￡斗￡ＣＰ咒ｆＰｒｏ，ＢＧＰ，ＣＮＰＣ，Ｚ＾“ｏｚ＾ｏ“０７２７５０，Ｃｈｉ靠口）

Ａｂｓｔｒ扯ｔＷａｖｅ

ｅｑｕａｔｉｏｎｆｏｒｗａｒｄｍｏｄｅｌｉｎｇｉｓ

ｏｎ

ａｖｅｒｙｉｍｐｏｒｔａｎｔｗａｙ

ｔｏ

ｓｔｕｄｙｔｈｅｓｅｉｓｍｉｃｐｒｏｐａｇａｔｉｏｎｕｎｄｅｒｇｒｏｕｎｄ．Ｉｎ

ｔｈｉｓｐａｐｅｒ，ｓｅｉｓｍｉｃｆｏｒｗａｒｄｍｏｄｅｌｉｎｇｂａｓｅｄｄｏｎｅ

ｕｓｉｎｇ

ｔｈｅ

ｓｔａｇｇｅｒｅｄ—ｇｒｉｄ

ｏｎ

ｔｗｏ－ｗａｙａｃｏｕｓｔｉｃｗａｖｅｅｑｕａｔｉｏｎａｎｄｔｗ０－ｗａｙｅｌａｓｔｉｃｗａＶｅｅｑｕａｔｉｏｎｉｓ

ｍｅｔｈｏｄ．ｍａｄｅ

Ｂｅｓｉｄｅｓ，

ａｃｃｏｒｄｉｎｇ

ｔ０

ｈｉｇｈ＿ｏｒｄｅｒ

ｄｉｆｆｅｒｅｎｃｅｔｈｅ

ｐｏｓｉｔｉｏｎｉｎｇｐｒｉｎｃｉｐｌｅ，

Ｉｎｔｈｅ

ｓｉｍｕｌａｔｉｏｎｂａｓｅｄ

ｐｒｏｇｒｅｓｓ

ｏｎｅ－ｗａｙａｃｏｕｓｔｉｃｗａｖｅｅｑｕａｔｉｏｎｉｓ

ｂｙｔｈｅＦｏｕｒｉｅｒｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅ

ｐｒｏｐａｇａｔｏｒ．

ｔｏ

ｏｆｃｏｍｐｕｔｉｎｇ，ｅａｃｈｇｅｏｐｈｏｎｅｄｉｒｅｃｔｌｙｒｅｃｅｉｖｅｓｔｈｅｗａｖｅｆｉｅｌｄｔｈａｔｈａｓｂｅｅｎｅｘｔｒａｐｏｌａｔｅｄ

Ｔｈｅｔｈｒｅｅｓｉｍｕｌａｔｉｏｎ

ｒｅｓｕｌｔｓｓｈｏｗ

ｔｈｅｓｕｒｆａｃｅｏｆｔｈｅ

ｇｅｎｅｒａｔｅ

ａ

ｍｏｄｅｌｗｉｔｈｏｕｔｌａｔｅｒａｌｓｔａｃｋｉｎｇ．ｗａｖｅａｎｄ

ｔｈａｔｔｗｏ．ｗａｙｗａ、ｒｅｅｑｕａｔｉｏｎｓ

ｄｉｒｅｃｔｉｖｅ

ｍｕＩｔｉｐｌｅｗａｖｅｗｈｉｃｈｒｅｄｕｃｅｔｈｅＳＮＲｓｅｖｅｒｅｌｙ，ｗｈｉｌｅｔｈｅｏｎｅ—ｗａｙｗａｖｅｅｑｕａｔｉｏｎｍｅｔｈｏｄｇｉｖｅｓ

ｃｒｅａｔｅｓｃｏｓｔｓ

ｈｉｇｈｅｒＳＮＲ

ｗｉｔｈｏｎｌｙｒｅｎｅｃｔｅｄｗａｙｅ，ｂｕｔｉｔｅｌａｓｔｉｃｗａｖｅｅｑｕａｔｉｏｎ

ｌａｒｇｅｒ

ｅｒｒｏｒｓ

ｗｉｔｈｗｉｄｅ－ａｎ一ｅｉｎｃｉｄｅｎｔｗａ’ｒｅ．Ｆｏｒｔｈｅｓａｍｅｍｏｄｅｌ，ｔｗｏ—ｗａｙ

ｃｏｓｔｓ

ａｌｇｏｒｉｔｈｍ

ｔｈｅ

ｍｏｓｔｃｏｍｐｕｔａｔｉｏｎａｌｔｉｍｅ，ｏｎｅ＿ｗａｙａｃｏｕｓｔｉｃｗａｖｅｅｑｕａｔｉｏｎａｌｇｏｒｉｔｈｍ

ａｌｇｏｒｉｔｈｍｉｎｂｅｔｗｅｅｎ．Ｓｏｆｏｒ

ａ

ｔｈｅｌｅａｓｔａｎｄｔｗ０－ｗａｙａｃｏｕｓｔｉｃｗａｖｅｅｑｕａｔｉｏｎｍｏｄｅｌｉｎｇＫｅｙｗｏｒｄｓ

ｃｏｍｐｌｅｘ

ｓｔｒｕｃｔｕｒｅｍｏｄｅｌ，

ｔｈｅｔｈｒｅｅ

ｍｅｔｈｏｄｓ

ｃａｎ

ｂｅｕｓｅｄｓｉｍｕｌｔａｎｅｏｕｓｌｙｗｉｔｈｄｉｆｆｅｒｅｎｔａｄｖａｎｔａｇｅｓａｎｄｄｉｓａｄｖａｎｔａｇｅｓ．

ｈｉｇｈ．ｏｒｄｅｒｄｉｆｆｅｒｅｎｃｅｍｅｔｈｏｄ，ｐｏｓｉｔｉｏｎｉｎｇｐｒｉｎｃｉｐｌｅ，

ｗａｖｅｅｑｕａｔｉｏｎｆｏｒｗａｒｄｍｏｄｅｌｉｎｇ，ｓｔａｇｇｅｒｅｄｇｒｉｄ，

ｏｎｅ．ｗａｙｗａｖｅｅｑｕａｔｉｏｎ

收稿日期作者简介

２００８＿０６一１２；

修回日期

２００８一０９－２ｚ．

皮红梅，１９８１年生，博士研究生，２００３年毕业于吉林大学地球物理专业，获理学学士学位，现在该校地球探测科学与技术学院攻读固体地球物理学专业博士学位，研究领域为复杂地震波场数值模拟．（Ｅ—ｍａｉｌ．ｂｒｏｏｋ＿３９９３＠ｅｍａｉｌ．ｊｌｕ．ｅｄｕ．ｃｎ）

通讯作者刘财，吉林大学教授，博士生导师，主要从事固体地球物理研究．（Ｅ－ｍａｍｌｉｕｃａＵｌｕ＠１６３．ｃｏｍ）

万方数据　

３９２

地球物理学进展

２４卷

Ｏ

引言

算效率低，正演记录的信噪比低．贺振华、熊高君（１９９８）等提出检波点下延原理［８￣儿］，首次采用单程声波方程ｍｄ３３将炮点和检波点的地震波场同时从地面向地下进行延拓得到共炮正演记录，但该方法是从偏移的思路来进行正演，计算效率较低，于是又提出一种能快速合成共炮记录的固定检波点定位原理ｎ引，此后，贺振华、胡光岷（２００４）等人提出的数学检波器的概念［１５］和贺振华、熊晓军（２００５）提出的等时叠加原理［１６１，从理论上对单程波动方程叠前正演方法进行了归纳和总结．本文在分析固定检波点定位原理的基础上，对其计算过程稍作修改，将延拓到地面的波场直接由每个检波点接收，无需横向叠加过程，得到了单程声波方程共炮记录，对三种波动方程得到的共炮记录进行分析对比．最后进一步分析了双程声波方程和单程声波方程的自激自收记录．

１

地震数值模拟方法是地震勘探和地震学的重要基础，是根据已知的地下介质结构和相应物理参数模拟研究地震波在地下各种介质中的传播规律，并计算在地面或地下各观测点的地震记录的一种地震模拟方法．数值模拟方法主要有三大类：波动方程法，积分方程法和射线追踪法．波动方程数值解法是建立在以弹性或粘弹性理论和牛顿力学为基础的双曲型偏微分方程一地震波传播方程的理论基础上的．由于地下介质性质不同，其相应的地震波传播方程也不同．如声学介质中的声波波动方程；弹性介质中的弹性波波动方程；粘弹性介质中的粘弹性波波动方程；孔隙弹性介质（双相或多相介质）中的双相或多相介质弹性波方程；各向异性介质中的各向异性弹性波波动方程等．积分方程法是建立在以惠更斯原理为基础的波叠加原理基础上的，其数学表达形式为波动方程的格林函数域积分方程式和边界积分方程式．射线追踪法是建立在以射线理论为基础的波动方程高频近似理论基础上的，其数学表达形式为程函方程和传输方程．

波动方程法又分为双程波方程法和单程波方程法．本文采用交错网格高阶有限差分法对双程声波方程和双程弹性波方程进行数值模拟．交错网格高阶有限差分法［１￣３］是董良国（２０００）首先提出的，他将Ｍａｄ耐ａｇａ（１９７６）提出的交错网格法［４３与Ｄａｂｌａｉｎ（１９８６）提出的高阶差分方法［５］有机结合，用于求解

双程声波方程数值模拟

直角坐标系下，各向同性介质中的二维声波方

程的一阶压力一速度方程ｍ３形式为：

警一一Ｋ（篆＋碧），

垫一一上丝．

ａ￡

ｐ

ａｚ’

婺一一土挈．

ａｔ

ｐ

＜１）

、１

ａｚ’

式中，户一户（ｚ，ｚ，￡）为声波波场，％＝砚（ｚ，２）为质点振动速度的水平分量，让一皱（ｚ，名）为垂直分量，Ｋ一即２为弹性模量，口为介质的纵波速度，ｌＤ为介质密度．

采用交错网格（如图１所示）高阶有限差分方法对方程（１）进行数值离散，其２阶时间２Ｎ阶空间的离散格式为：

横向各向同性介质一阶速度一应力弹性波方程，其

后他又对该方法的稳定性问题［６］（２０００）和频散问题＂］（２００４）进行了研究．双程波方程法包含了地震波场的全部信息，可以处理复杂的地质构造，但是计

１．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．一——．．——————一．．——．——．．

矿％∽＝一∽沪警｛塾Ⅲ＋竽，Ｊ）一谚（ｚ一竿，歹）］

十塾［破（奶＋学）一破（叫一竽）］｝，

谚

ｌ

＋

１—２

一记

．Ｚ

＋

１—２

一

妻“［户蚪｝（ｉ＋咒，歹）一户外专（ｉ＋１一行，歹）］，（２）

矿

ｌ

｝

１—２一

破

．￡

｝

１—２一

其中：∥＋专（ｉ，＿『）表示户（ｚ，ｚ，￡＋譬），

万方数据　

塑一

Ⅳ

惫箭一∞一∞

∑ｃ。［户蚪｝（ｉ，歹＋咒）一夕蚪｝（ｚ，歹＋１一咒）］．

＾＝１

谚（ｆ＋丢，＿ｆ）表示Ｖ。（ｚ＋等，２，ｔ），其它与此类

２期皮红梅，等：波动方程数值模拟的三种方法及对比

３９３

同．厶为交错网丰子差分系识．

２双程弹性波方程数值模拟

均匀各向同性介质中的二维弹性波方程ｍ３为：

．ｏ等＝鲁＋等，Ｐ警＝鲁＋等，

鲁＝（Ａ＋２ｐ）卺＋｜；Ｉ卺，

鲁＿＿（Ａ＋刎凳＋Ａ警，鲁＝ｐ（等＋卺）．

其中，让和钞。分别为质点振动速度的水平分量和垂直分量；％，％，屹为应力分量；ｌＤ为介质密度；Ａ和口为拉梅常数．

■ｐｘ

１，

●ｖｘ’，ｐ

∥

ｏ屹，形

图１声波方程交错网格示意图

Ｆｉｇ．１

Ｔｈｅｓｔａｇｇｅｒｅｄｇｒｉｄｄａｃｏｕｓｔｉｃｗａｖｅｅｑｕａｔｉｏ矗

采用交错网格（如图２所示，据Ｌｅｖａｎｄｅｒ［１９３）高阶有限差分法来求解方程（３）．

图２

弹性波方程交错网格示意图（据Ｌｅｖａｎｄｅｒ１９８８）

Ｆ螗．Ｚ

Ｔｈｅｓｔａｇｇｅｒｅｄｇｒｉｄｏｆｅｌａｓｔｉｃ

ｗａｖｅｅｑｕａｔｉｏｎ（Ｌｅｖａｎｄｅｒ１９８８）

万　

方数据呓｛＝吒吾＋考惫＊［蚤“《字。，一Ｐ：２孚≯

＋∑ｃ。（《．升学一≤，卜学）］，

’（４）

咙Ｕ２，尸：。／２．』，ｓ：一。，。分别为质点速度分量饥和应

力分量盯。，％的离散值，其它几个分量的离散方法与此类似．

３单程声波方程数值模拟

二维各向同性均匀介质中的单程波方程田３为：

筹一土ｉ志庐＝士ｉ√等一磋夕，，

（５）

歹＝声（走。，名，甜）为频率一波数域声波波场，“＋”、“一”号分别代表上行波方程和下行波方程．则根据Ｄ．Ｒｉｓｔｏｗ（１９９４）提出的傅立叶有限差分‘２１３（ＦＦＤ）法可得用于正演的波场延拓算子为：

ｆ声（忌。，锰ｌ，ｃ￡，）一声（惫：，≈，（ｃ，）ｅ＿‘√孟巧△ｆ

｛石（ｚ，铣ｌ，叫）一万（ｚ，巧，（ｃ，）ｆｂ血ｂ∞缸

．

（６）

【户（ｚ，ｚ灶１，（ｃ’）＝夕（ｚ，≈，ｃＥ，）ｒ巩缸

其中：石＝石（ｚ，２，∞）为频率一空间域声波波场，ｃ为

每一个延拓深度层内的背景速度，取为：

ｃ＝ｍｉｎ（口（ｚ，２。）），

纵…，一志一÷，

Ａ＿钟吲【蘸Ｊ’

ｆ生竺］

６＝丢［（詈）２＋詈＋－］．

根据贺振华，熊高君提出的固定检波点定位原理［Ｈ３合成共炮记录，对其计算过程稍作修改，将延拓到地面的波场直接由每个检波点接收，无需横向叠加过程．原理表述如下：固定炮点和检波点，用下行波方程将震源波场从炮点位置延拓到反射界面，与所有反射点的位置做相关运算，相当于模拟入射波沿反射界面激发新震源或反射波，把这个新震源暂时置于所产生的位置，同时把原来的入射波继续向下延拓直到最深处，形成所有深度处的新震源，再从最深处用上行波方程将新震源向上延拓，每延拓一步，将延拓波场与该处的新震源波场叠加，将叠加后的波场继续向上延拓，直到检波器所在位置，将该延拓波场与检波器谐振因子进行相关，该相关结果就是所固定检波点的正演结果．

记：地面记录５９（ｚ，￡），震源波场ｓ（ｚ，２，ｆ），检

其２阶时间、２Ｎ阶空间的离散格式为：

地球物理学进展

２４卷

波谐振波场ｇ（ｚ。，ｚｏ，ｔ），反射系数ｒ（ｚ，ｚ）．其中：

ｓ（ｚ，州）一Ｆ五嘞一几ｚ一五

（７）

ｌＯ，ｚ≠ｚ，

（ｚ．，勐）、（ｚ。，‰）分别为震源位置和检波点位置，本论文中震源和检波因子都用ｒｉｃｋｅｒ子波．将震源和检波因子都变换到频率域ｉ（ｚ，ｚ，ｃｃ，）、虿（以，孙，甜），并且反射系数拓展为频率不变函数尸（ｚ，名，∞），即

尹（ｚ，ｚ，∞）＝ｒ（ｚ，２）．

固定检波点原理Ｐ－Ｐ波合成记录的步骤如下：（１）根据（６）式将震源波场向下延拓，直到地下最深处，并且与延拓深度上的反射系数进行相关，形成新震源：

ｊ（ｚ，ｚ件１，∞）＝ｊ（ｚ，ｚ，，（￡’）ｅ一让。血，

（８）西（ｚ，≈＋ｌ，叫）一ｊ（ｚ，乃＋１，（￡，）尹（ｚ，影件１，叫）．

（９）

（２）根据（６）式从地下最深处将反射界面上的所有新震源向上延拓，并与该延拓层新震源叠加，直到将波场延拓到检波器所在的深处（假设地面ｚ＝０处）．

歹１（ｚ，争１，ｃ￡，）＝歹ｌ（ｚ，乃，倒）ｅ－让ｒ缸＋而（ｚ，夸ｌ，（ｃＪ）

（１０）

其中歹ｌ（ｚ，名。，叫）一鄢（ｚ，名。，ｃ￡，）

（３）将延拓到地面的波场歹，（ｚ，ｚ＝Ｏ，∞）与检波器对应位置的值乘以频率域检波因子，即为该检波器接收到的波场值：

元（ｚｇ，∞）＝歹１（ｚＦ，ｚ—Ｏ，ｃｃ’）虿（ｚｇ，ｚ，∞），

（１１）

（４）将合成记录变换到时间域，得到以（ｚ，￡）．固定新的检波点位置，重复步骤（３）和（４），可得到完整的共炮记录户（ｚ，ｆ）．

对于转换波Ｐ—ＳＶ波的合成，只需将步骤（２）中向上延拓的波速用横波速度代替即可．而对于零偏移距正演记录，根据爆炸反射界面原理［２引，只需将步骤（１）中的新震源用地下反射系数代替，再进行步骤（２）和（３）的运算，即可获得．

４数值计算结果

图３为逆断层网格模型，介质泊松比为Ｏ．２５，模型大小为５０００

ｍ×３０００

ｍ，离散网格间距为如

一如＝１０ｍ，震源子波为主频２５Ｈｚ的雷克子波，位于地表２５００ｍ处，采用的交错网格高阶有限差分法精度为时间２阶空间８阶，分别进行双程声波波动方程模拟和弹性波方程模拟，合成共炮记录如图４和图５所示．单程波模拟的时间采样率为ｌｍｓ，记录长度为２．５ｓ，Ｐ—Ｐ波合成记录和Ｐ—ＳＶ波合成记录如图６所示．对比三种情况下的合成记录

万　

方数据可以看出，弹性波波场成份丰富，较正确地反映了波在界面上的反射、透射等物理现象，符合地震波的传

播规律，图５中Ｓ波双曲线顶部出现较明显的频散现象，这主要是由于弹性波方程采用了与声波方程相同的离散参数所致，所以波场模拟需要考虑Ｓ波时，应当以Ｓ波频散最小为原则．图６单程波模拟结果不含直达波、多次波等干扰波，只含有与地质界面相对应的一次反射波，对比图５、６双曲线两翼可以发现，单程波波场能量主要集中在小炮检距范围内，这是因为单程波方程是双程波方程的近似，单程波算子存在一个最大倾角，对该倾角范围之外的波模拟误差较大．

图４双程声波共炮记录

Ｆｉｇ．４

Ｓｈｏｔｒｅｃｏｒｄｓｏｆｔｗｏ—ｗａｙａｃｏｕｓｔｉｃｗａＶｅｅｑｕａｔｉｏｎ

图７为地层尖灭网格模型，模型大小为５０００

ｍ

×３０００

ｍ，离散网格间距为ｄｚ＝如＝５ｍ，时间采样

率为１ｍｓ，记录长度为２．２ｓ，震源子波为主频２５Ｈｚ的雷克子波，双程声波方程和单程波方程自激自

收记录如图８所示．对比图８两种记录可见，单程波

２期皮红梅，等：波动方程数值模拟的ｊ种方法及对比

３９５

方程记录各反射界面的反射波很明显，而双程声波受起伏界面影响，双程声波模拟产生的多次波将倾方程记录倾斜地层反射波模糊不清，这主要是由于

斜地层反射波掩盖所致．

图７地层尖灭网格模型

Ｆｉｇ．７

Ｓｔｒａｔａｐｉｎｃｈｉｎｇ

ｍｏｄｅｌ

万　

方数据

３９６

地球物理学进展

２４卷

图８自激自收记录

左：双程声波方程，右：单程声波方程

Ｆｉｇ．８

ｌｅｆｔ：ｔｗｏ．ｗａｙ

ｗａｖｅ

Ｚｅｒｏｏ“ｓｅｔｒｅｃｏｒｄｓ

ｗａＶｅｅｑｕａｔｉｏｎ

ｅｑｕａｔｉｏｎ，ｒｉｇｈｔ：ｏｎｅ－ｗａｙ

Ｗａｎｇ

Ｘ

Ｍ，ＺｈａｎｇＨ

ｉｎ

Ｌ，ＷａｎｇＤ．Ｍｏｄｅｌｉｎｇ

ｏｆ

ｓ幽ｍｉｃ

ｕｇｉｎｇ

ａ

ｗａ、ｒｅ

５

结论

ｐｒｏｐａｇａｔｉｏｎ

ｏｒｄｅｒ

ｈｅｔｅｍｇｅｎｅｏｕｓ

ｐｏｒｏｅｌａｓｔｉｃｍｅｄｉａ

ｈｉｇｈ－Ｊ．

基于不同波动方程数值模拟结果分析表明，双程波方程数值模拟的波场成份丰富，比较实际地反映了地震波在介质中传播的基本情况，但是当介质比较复杂时，多次波严重干扰了界面的一次反射波信息，波型识别困难．单程声波数值模拟结果可以清晰地反映地质分界面，但是受倾角限制的影响，波场能量主要集中在小炮检距范围内．因此对于复杂地质模型，三种模拟方法可以取长补短，综合应用．致谢感谢东方地球物理公司物探技术研究中心方法部田振平，卢秀丽等同志在论文完成过程中给予的热情帮助．参

考

文

［４］

ｓｔａ昭ｅｒｅｄｆｉｎｉｔＰｄｉｆｆｅｒｅｎｃｅｍｅｔｈｏｄ［Ｊ］．

Ｃｌｌｉｎｅｓｅ

Ｇｅｏｐｈｙｓ．（ｉｎＣｈｉｎｅｓｅ），２００３，４６（６），８４Ｚ～８４９．ＭａｄａｒｉａｇａＲ．

ＤｙｎａｒＩｌｉｃｓｏｆ

ａｎ

ｅｘｐａｎｄｉｎｇｃｉｒｃｕＩａｒｆａｕｌｔ［Ｊ］．

ＢＳＳＡ，１９７６，６６（３）：６３９～６６６．

［５］ＤａｂｌａｉｎＭ

ｓｃａｌａｒｗａｖｅ

Ａ．Ｔｈｅ

ａｐｐｌｉｃａｔｉｏｎｏｆ

ｈｉｇｈ—ｏｒｄｅｒ

ｄｉｆ｛ｅｒｅｎｃｉｎｇ

ｔｏ

ｔｈｅ

ｅｑｕａｔｉｏｎ［Ｊ］．Ｇｅｏｐｈｙｓｉｃｓ，１９８６，５１（１）１５４～６６．

［６］董良国，马在田，曹景忠．一阶弹性波方程交错网格高阶差分解

法稳定性研究口］．．地球物理学报，２０００，４３（６）：８５６～８６４．

Ｄｏｎｇ

Ｌ

Ｇ。ＭａＺＴ，ＣａｏＪＺ．Ａ

ｈｉｇ卜ｏｒｄｅｒ

ｓｔｕｄｙ

ｏｎ

ｓｔａｂｉｌｉｔｙｏｆｔｈｅｏｆ

ｆｉｒｓｔ－ｏｒｄｅｒ

ｓｔａｇｇｅｒｅｄ—ｇｒｉｄｅＩａｓｔｉｃｗａｖｅ

ｄｉｆｆｅｒｅｎｃｅｍｅｔｈｏｄ

ｅｑｕａｔｉｏｎ［Ｊ］．ｃｈｉｎｅｓｅＪ．Ｇｅｏｐｈｙ３．（ｉｎｃｌｄｎｅｓｅ），

ＺＯＯＯ。４３（６）：８５６～８６４．

［７］董良国，李培明．地震波传播数值模拟中的频散问题［Ｊ］．天然

气工业，２００４，２４（６）：５３～５６．

Ｄ０ｎｇ

Ｌ

Ｇ，ＬｉＰＭ．

Ｄｉｓｐｅｒｓｉｖｅ

ｐｒｏｂｌｅｍｉｎ

ｓｅｉ８ｍｉｃ

ＷａＶｅ

献（Ｒｅｆｅｒｅｎｃｅｓ）：

［８］

ｐｒｏｐａｇａｔｉｏｎｎｕｍｅｒｉｃａｌ２００４，２４（６）：５３～５６．

ｒＩｌｏｄｅ“ｎｇ［Ｊ］．ＮａｔｕｒａｌＧａｓ

Ｉｎｄｕｓｔｒｙ，

［１］董良国，马在田，曹景忠，王华忠，耿建华，雷兵，许世勇．一阶弹

性波方程交错网格高阶差分解法［Ｊ］．地球物理学报，２０００，４３

（３）：４１１～４１９．Ｄｏｎｇ

Ｌ

Ｈｅ

ＺＨ，ＸｉｏｎｇＧ

Ｍｏｄｅｌｉｎｇ

Ｊ，ｚｈｎｇＹＸ．Ｎｏｌｌｚｅｒｏｏ“ｓｅｔＳｅｉｓｍｉｃＢｙ∞ｅ．ｗａｙＡｃｏｕ“ｃｗａｖｅ－ｅｑｕａｔｉｏｎ［Ｊ］．

Ｆｏ蝴８ｒｄ

ＳＥＧＥｘｐａｎｄｅｄＡｂｓｔｒａｃｔｓ，１９９８．

Ｇ，ＭａＺＴ，ＣａｏＪＺ，以ｎＺ．Ａｓｔａｇｇｅｒｅｄ－ｇｒｉｄｈｉｇｈ－

ｗａｖｅ

［９］熊高君，贺振华，张琳，等．共炮记录正演模拟检波点下延记录

原理［Ｊ］．石油物探，１９９９，３８（２）：４３～４９．

ＸｉｏｎｇＧ

ｏｒｄｅｒｄｉｆｆｅｒｅｎｃｅｍｅｔｈｏｄｏｆｏｎｅ－ｏｒｄｅｒｅｌａｓｔｉｃ

ｅｑｕａｔｉｏｎ［刀．

ＣｈｉｎｅｓｅＪ．Ｇｅｏｐｈｙｓ．（ｉｎＣｈｉｎｅｓｅ），２０００，４３（３）：４１１～４１９．

Ｊ，ＨｅＺＨ，ＺｈａｎｇＬ。甜口Ｚ．Ｇｅｏｐｈｏｎｅｒｅｃｏｒｄｆｒｏｍ

ｃｏｍｍｏｎ－８ｈｏｔ

”ｃｏｒｄ

ｆｏ聊ａｒｄ

［２］汪功礼，张庚骥，崔锋修，等．三维感应测井响应计算的交错罔

格有限差分法［Ｊ］．地球物理学报，ｚ００３，４６（４），５６１～５６７．Ｗａｎｇ

Ｇ

ｍｏｄｅｌｉｎｇ

［Ｊ］．ＧｅｏｐｈｙｓｉｃａＩ

ＰｒｏｓｐｅｃｔｉｎｇｆｏｒＰｅｔｒｏｌｅｕｍ，１９９９，３８（２）１４３～４９．

Ｌ，Ｚｈａｎｇ

ｇ—ｄ

ＧＪ，ＣｕｉＦ

Ｘ，ｃｔ“．Ａｐｐｌｉｃａｔｉｏｎｏｆ

ｔｏ

［１０］熊高君，贺振华，黄德济．等．转换横波共炮记录单程声波方程

正演模拟［ｊ］．物探化探计算技术，１９９８，２０（３）：２１３～２１７．

Ｘｉｏｎｇ

ｓｔａｇｇｅｒｅｄ

３一Ｄ

ｆｉｌｌｉｔｅｄｉ“ｅｒｅｎｃｅｍｅｔｈｏｄ

ｔｈｅｃｏｍｐｕｔａｔｉｏｎｏｆ

ｉｎｄｕｃｔｉｏｎ１０９９ｉｎｇｒｅｓｐｏｎｓｅ［Ｊ］．ＣｈｉｎｅｓｅＪ．Ｇｅｏｐｈｙｓ．（ｉｎＧＪ，ＨｅＺ

Ｈ，Ｈｕａｎｇ

ｏｆ

Ｄ

Ｊ，“以．（：ｏｍｍｏｎ—ｓｈｏｔｒｅｃｏｒｄｓ

ｏｎｅ＿ｗａｙ、张Ｖｅ

Ｃｈｉｎｅｓｅ），２００３，４６（４）：５６１～５６７．

ｆｏｒｗａｒｄｍｏｄｅｌｉｎｇ

ｃｏｎｖｅｒｔｅｄｓｈｅａ卜ｗａＶｅｗｉｔｈ

ｔｅｃｈｎｉｑｕｅｓ

［３］

王秀明。张海澜，王东．利用高阶交错网格有限差分法模拟地震波在非均匀孔隙介质中的传播［Ｊ］．地球物理学报，２００３，４６

（６），８４２～８４９．

ｅｑｕａｔｉｏｎ［Ｊ］．（ｂｍｐｕｔｉｆｌｇ

ｆｏｒ

ｇｅｏｐｈｙｓｉｃａｌａｎｄ

ｇｅｏｃｈｅｍｉｃａｌｅｘｐｌｏｒａｔｉｏｎ，１９９８，２０（３）：２１３～２１７．

［１１］熊高君，贺振华，黄德济ｔ等．多次反射可控制的波动方程共

万方数据　

２期皮红梅，等：波动方程数值模拟的三种方法及对比

３９７

炮记录正演模拟口］．石油地球物理勘探，１９９８，３３（５）：５８３～

５９０．

Ｘｉｏｎｇ

Ｇ

Ｊ，ＨｅＺＨ，Ｈｕａｎｇ

Ｄ

Ｊ，以以．ＷａＶｅ－ｅｑｕａｔｉｏｎ

ｃｏｍｍｏｎ－ｓｈｏｔｓｅｉｓｍｏｇｒａｍｆｏｒｗａｒｄｓｉｍｕｌａｔｉｏｎ

ｗｉｔｈ

ｍｕｌｔ卜

ｒｅｆｌｅｃｔｉｏｎ

ｗａｖｅ

ｃｏｎｓｔｒａｉｎｅｄ［Ｊ］．０ｉｌＧｅｏｐｈｙｓｉｃａｌ

Ｐｒ０８ｐｅｃｔｉｎｇ，１９９８。３３（５）：５８３～５９０．

［１２］刘札农，崔风林，张剑锋．三维复杂构造中地震波模拟的单程

波方法［Ｊ］．地球物理学报，２００４，４７（３）；５１４～５２０．

Ｌｉｕ

Ｌ

Ｎ．ＣｕｉＦ

Ｌ。ＺｈａｎｇＪＦ．Ｓｅｉｓｍｉｃ

ｎｌｏｄｅｌｉｎｇ

ｗｉｔｈ

ｏｎｅ－

ｗａｙ、张ｖｅ

ｅｑｕａｔｉｏｎｉｎ

３Ｄ

ｃｏｍｐｌｅｘ

ｓｔｒｕｃｔｕｒｅｓ．ＣｈｉｎｅｓｅＪ．Ｇｅｏｐｈｙｓ．（ｉｎＣｈｉｎ舶ｅ）ｔ２００４，４７（３）ｉ５１４～５２０．

［１３］谢桂生，刘洪。李幼铭，等．界面起伏条件下反射／透射算子＋

单程波方程的地震波模拟方法［Ｊ］．地球物理学报，２００５．４８

（５）：１１７２～１１７８．Ｘｉｅ

Ｇ

Ｓ．“ｕＨ，Ｌｉ

Ｙ

Ｍ，以以．Ｓｅｉｓｍｉｃｍｏｄｅｌｉｎｇ

ｂｙｒｅｎｅｃｔｉｏｎ／ｔｒａｎｓｍｉｓｓｉｏｎｏｐｅｒａｔｏｒ

ａｎｄ

ｏｎｅ－ｗａｙ

ｗａｖｅ

ｅｑｕａｔｉｏｎ

ｕｎｄｅｒｔｈｅｃｏｎｄｉｔｉｏｎ

ｏｆ

ｆｌｕｃｔｕａｔｉｎｇ

ｒｅｆｌｅｃｔｏｒｓ［Ｊ］．ＣｈｉｎｅｓｅＪ．

Ｇｅｏｐｈｙｓ．（ｉｎＣｈｉｎｅｓｅ），２００５．４８（５）：１１７２～１１７８．

［１４］熊高君，贺振华，黄德济，等．改进的正演模拟定位原理［Ｊ］．

石油地球物理勘探，１９９８，３３（６）：７４２～７４８．

ＸｉｏｎｇＧＪ，ＨｅＺＨ，ＨｕａｎｇＤ

Ｊ，以“．１ｍｐｒｏｖｅｄ

ｐｏｓｉｔｉｏｎｉｎｇ

ｐｒｉｎｃｉｐｌｅ

ｏｆ

ｆｏ唧ａｒｄ

ｍｏｄｅＩｉｎｇ［Ｊ］．ｏｉｌ

ＧｅｏｐｈｙｓｉｃａｌＰｒｏｓｐｅｃｔｉｎｇ，１９９８，３３（６）：７４２～７４８．［１５］贺振华。胡光岷，黄德济．数学检波器与波动方程地震叠前正演［Ｊ］．成都理工大学学报（自然科学版），ｚ００４，３ｌ（６）：６７５～６７８．

ＨｅＺＨ，ＨｕＧＭ，Ｈｕａｎｇ

Ｄ

Ｊ．Ｍａｔｈｅｍａｔｉｃａｌｇｅｏｐｈｏｎｅｓｗｉｔｈａｐｐｌｉｃａｔｉｏｎｓｔｏ

ｓｅｉｓＩＴｌｉｃ

ｐｒｅｓｔａｃｋｉｎｇｆｏｒｗａｒｄｍｏｄｅｌｉｎｇ［Ｊ］．

ＪｏｕＭｎｌｏｆ

ＣｈｅｎｇｄｕＵｎｉｖｅｒｓｉｔｙ

ｏｆＴｅｃｈｎｏｌｏｇｙ（Ｓｃｉｅｎｃｅ＆

万　

方数据ＴｅｃｈｎｏｌｏｇｙＥｄｉｔｉｏｎ），２００４，３１（６）：６７５～６７８．

［１６］贺振华，熊晓军．等时叠加波动方程叠前正演［Ｊ］．物探化探

计算技术，２００５，２７（３）：１９４～２００．

ＨｅＺＨ，ＸｉｏｎｇＸＪ．Ｓｅｉｓｎｌｉｃｐｒｅｓｔａｃｋｉｎｇｆｏ邢ａｒｄ

ｍｏｄｅｌｉｎｇ

ｂａｓｅｄ

ｏｎ

ｅｑｕａｌ－ｔｉｍｅ

ｓｔａｃｋ

ｐｒｉｎｃｉｐｌｅ

ａｎｄ

ｏｎｅ－ｗａｙ

ｗａｖｅ

ｅｑｕａｔｉｏｎ［Ｊ］．ｃｏｍｐｕｔｉｎｇＴｅｃｈｎｉｑｕｅｓｆｏｒＧｅｏｐｈｙｓｉｃａｌａｎｄ

ＧｅｏｃｈｅｍｉｃａｌＥｘｐｌｏｒａｔｉｏｎ，２００５，２７（３）：１９４～２００．

［１７］牟永光，裴正林著．三维复杂介质地震数值模拟［Ｍ］．北京：

石油工业出版社，２００４．

Ｍｕ

Ｙ

Ｇ，ＰｅｉＺＬ．Ｓｅｉｓｍｉｃｎｕｍｅｒｉｃａｌｍｏｄｅｌｉｎｇ

ｆｏｒ３－ｄ

ｃｏｍｐｌｅｘ

ｍｅｄｉａ［Ｍ］．ＢｅｉｊｉｎｇｆＰｅｔｒｏｌｅｕｍＩｎｄｕｓｔｒｙＰｒ船ｓ，２００４．

［１８］Ｖ试ｅｕｘＪ．Ｐ—ｓＶｗａｖｅ

ｐｒｏｐａｇａｔｉｏｎ

ｉｎ

ｈｅｔｅｒｏｇｅｎｅｏｕｓｍｅｄｉａ：

Ｖｅｌｏｃｉｔｙ＇ｓｔｒｅｓｓｆｉｎｉｔｅ－ｄｉｆｆｅｒｅｎｃｅｍｅｔｈｏｄ［Ｊ］．

Ｇｅｏｐｈｙｓｉｃｓ。

１９８６，５１（４）；８８９～９０１．

［１９］Ｌ＿ｅｖａｎｄｅｒ

Ａ

Ｒ．

Ｆｏｕｒｔｂｒｏｒｄｅｒｆｉｎｉｔｅ＿ｄｉｆｆｅｒｅｎｃｅＩＬｓＶ

ｓｅｉｓｍｏｇｒａｍｓ［Ｊ］．Ｇｅｏｐｈｙｓｉｃｓ，１９８８，５３（１１）：１４２５～１４３６．

［２０］马在田著．地震成像技术一有限差分法偏移［Ｍ］．石油工业出

版社．１９８９．

ＭａＺ

Ｔ．

Ｓｅｉｓｍｉｃｉｍａｇｉｎｇｔｅｃｈｎｉｑｕｅｓ－Ｍｉｇｒａｔｉｏｎ

ｗｉｔｈｆｉｎｉｔｅ

ｄｉｆｆｅｒｅｎｃｅ

ｍｅｔｈｏｄ［Ｍ］．Ｂｅｉｊｉｎｇ：ＰｅｔｒｏｌｅｕｍＩｎｄｕｓｔｒｙ

Ｐｒｅｓｓ，

１９８９．

［２１］

ＲｉｓｔｏｗＤ，ＲｕｈｌＴ．Ｆｏｕｒｉｅｒ“ｎｉｔｅ

ｄｉｆｆｅｒｅｎｃｅ

ｍｉｇｒａｔｉｏｎ［Ｊ］．

Ｇｅｏｐｈｙｓｉｃｓ，１９９４，５９（１２）：１８８２～１８９３．

［ｚ２］贺振华，王才经，李建朝，等．反射地震资料偏移处理与反演

方法［Ｍ］．重庆大学出版社．１９８９．

Ｈｅ

Ｚ

Ｈ，Ｗａｎｇ

ＣＪ，“Ｊ

Ｃ，甜“．Ｍｉｇｒａｔｉｏｎａｎｄ

ｃｏｎＶｅｒｓｉｏｎ

ｍｅｔｈｏｄｓｆｏｒ

ｒｅｆｌｅｃｔｉｏｎｓｅｉｓｍｉｃ

ｄａｔａ［Ｍ］．ｃｈｏｎｇｑｉｎｇ

Ｕｎｉｖｅｒｓｉｔｙ

Ｐｒｅｓｓ，１９８９．