波动方程的变步长有限差分数值模拟
收稿日期：２００７－０３－２３；修订日期：２００７－０４－２７
作者简介：李胜军，男，在读硕士研究生，研究方向为地震波传播理论。联系电话：（０５４６）８３９２０５５，Ｅ－ｍａｉｌ：ｈｄｐｕｌｉｓ＠１２６．ｃｏｍ，通讯地址：（２５７０６１）中国石油大学（华东）地球资信与信息学院。
＊中国石油大学（华东）研究生创新基金资助，编号：Ｓ２００６—０６。
油气地球物理
２００７年７月
ＰＥＴＲＯＬＥＵＭＧＥＯＰＨＹＳＩＣＳ
第５卷第３期
在地震资料采集、处理和解释中通常需要进行地震波场数值模拟：假设已知地下的地质情况，应用地震波运动学和动力学的基本原理，计算给定地质模型的地震响应。这种做法对正确认识地震波的运动学和动力学特征，以及准确分析油气藏的反射波场特征有着重要的指导意义。声波在介质中的正演模拟研究为我们精确模拟地震波在复杂介质中的传播提供了理论基础［１］。
傅立叶变换法和高阶有限差分法（ＦＤ）已成为计算声波方程空间导数的标准技术［２，３］。虽然常网格步长差分算法比较容易实现，但是它们对大部分模型都增大了不必要的计算量。例如，对存在浅层低速带的沉积盆地模型地面地震记录进行模拟时，由于低速地层阻抗小，地震波传入其中会引起较大的振幅和较长的延续时间（这与深层的高速层完全不同）。由于这些浅层低速层中地震波的波长较短、地层厚度较小，模拟时需要用小网格进行。这样，常网格步长算法就必须用小网格离散整个模型，从而增加了不必要的代价，如内存、计算量的增大。
因而，采用变网格算法将能改进有上覆低速层情况模拟结果的有效性（对地层中间有超薄夹层的情形，必须用精细网格覆盖才能精确的对地层进行模拟）。应用这种变网格算法既能实现对夹层的模拟，又能保障计算量不增加。因此这种通过函数实现在任意深度上网格步长变化的有限差分方法被
推广［４］。为了计算空间导数，在Ｘ方向用傅立叶变换法或有限差分算法，在Ｚ方向使用高阶有限差分方法。通过时间积分快速展开法（ＲＥＭ）来保障差分方法的计算精度［６］。这种差分技巧比二阶时间差分有较高的精确度且计算用时短。
１时间积分
均匀介质中的二维声波方程可用下式表示［２］
式中：Ｐ＝Ｐ（ｘ，ｚ，ｔ），代表压力项；ｃ＝ｃ（ｘ，ｚ），代表速度；ｓ＝ｓ（ｓ，ｚ），代表震源函数；Ｌ２为差分算子。在密度!＝!（ｘ，ｚ）变化的情况下，常用的是Ｖｉｄａｌｅ给出的公式［５］
波动方程的变步长有限差分数值模拟＊
李胜军１，２）
孙成禹１）
张玉华１）
倪长宽１）
１）中国石油大学地球资源与信息学院；２）中石油勘探开发研究院西北分院
摘要：有限差分算法是常用的正演模拟方法之一，其包含的地震信息丰富，且实现简单。传统的有限差分方法通常都采用均匀网格步长，在对含低速／高速介质、
薄层／厚层介质的模型进行波场模拟时往往缺乏稳定性。文章介绍了一种可以有效解决上述问题的变网格算法，对常规有限差分法与变网格差分算法在内存需求、计算速率等方面的差别进行了比较，对变网格差分算法中的边界条件、
时间积分的快速展开算法作了阐述，进而总结了变网格算法的优点。关键词：变步长；边界条件；计算时间；快速展开法；数值模拟
!２
Ｐ!ｔ２＝－Ｌ２Ｐ＋ｓ
（１）
（２）
－Ｌ２
＝ｃ
２
!２!ｘ２＋!２
!ｚ
２"
#
（３）
（４）
!２
Ｐ!ｔ
２＝－Ｌ２Ｐ!"$
－１!Ｌ２Ｐ＋ＰＬ２１!＋ｓ
－Ｌ２
＝!ｃ
２
２
!２!ｘ２＋!２
!ｚ
２%
$
波动方程的时间积分可通过Ｋｏｓｌｏｆｆ等提出的快速展开法（ＲＥＭ）来计算［６］。为简便起见，仅对基本方程（１）的快速展开法做一些介绍。对变密度的情况可用类似的方法进行推导。
对震源项ｓ（ｘ，ｚ，ｔ）＝ｇ（ｘ，ｚ）ｈ（ｔ），边界条件为吸收边界时，式（１）的常规解为
该常规解通过契比雪夫展开近似为
式中：!为变量；是空间差分算子Ｌ２的数值逼近；Ｑｋ指修改了的契比雪夫多项式；系数ｂｋ为
式中：Ｊｋ为ｋ阶精度的贝赛尔函数。
如果空间参数－的特征值很靠近实轴、Ｒ２比
最大特征值大，当Ｎ＞Ｒ?ｔ时，（６）、（７）式按指数规律收敛。参数－的最大特征值由Ｘ方向上的差分参数的最大特征值Ａｘ和Ｚ方向上的相应特征值Ａｚ的和来决定。这里，Ｒ２可由（８）式计算
２空间导数计算
通常情况下，勘探区域的地表较为复杂且介质中波的传播速度较低。为较好地压制数值频散并提高模拟结果的分辨率，就要在整个计算区域采用小网格步长来做差分数值模拟。这样就造成了对高速区域的过采样，浪费了计算机资源。为减小内存需求和节省计算时间，可以用图１所示方法进行模拟。对于地层中有低速层的情形可以采用图２所示方法来达到提高分辨率并且计算量增加不大的差分算法。
图３给出了简单的模型：用固定小网格步长!ｘ和"ｚ采样，可以很好地满足低速层的采样需求，但是这样造成了对高速层的过采样。为避免这种过采样，可以从某一深度ｚ０（图１）开始采用双倍的步长来采样。在修改的网格上，Ｘ方向的导数用傅立叶法或有限差分法计算，#ｘ通过采样定理来确定。由于
$ｘ在Ｘ方向上是连续变化的，用常规的有限差分
算法就可计算。但是在Ｚ方向上，%ｚ是不连续的，基
于连续变化的傅立叶变换法不再适用，所以Ｚ方向上的导数要用高阶有限差分法近似。
为了获得一个稳定的算法，使用对称的有限差分参数。记Ｐ（ｉ!ｘ，ｊ!ｚ）＝Ｐｉ，ｊ，一种用对称有限差分参数近似的二阶导数定义如下［７］
对应的双步长网格定义如下
式中：２Ｎ是算子长度［８］。
（５）
Ｐｘ，ｚ，!"ｔ＝
ｔ０
#ｓｉｎ!!
"ＬＬ
ｈｔ－!"!ｄ$%!ｇｘ，!"ｚ（６）
Ｐｘ，ｚ，!"ｔＮ／２
ｋ＝０&ｂ２ｋ＋１Ｒ
ｉＬ
!Ｑ２ｋ＋１ｉＬ
!Ｒ!"$
%
ｇｘ，
!"ｚＬ
!２（７）
ｂｋ＝
ｔ
０#
Ｊ２ｋ＋１!!"
ＲＲ
ｈｔ－!"!ｄ!Ｌ!２Ｌ
!２（８）
ｘ
ｚ
ｚ０
图１低速地表模型差分网格分布示意图
ｘ
ｚ
图２低速夹层模型差分网格分布示意图
ｚ
ｘ
１５００ｍ／ｓ
３０００ｍ／ｓ
图３二维模型
（９）
Ｒ２＝ｃ
２
Ａｘ
!!"ｘ２
＋Ａ
ｚ
!!"
ｚ２
!"ｍａｘ
（１０）
Ｄｚ２
$%Ｐｉ，ｊ＝"２
Ｐ"ｚ
２
ｉ，ｊ
＝１!!"ｚ２Ｎ
ｋ＝１&#ｋＰｉ，ｊ＋ｋ－２Ｐｉ，ｊ＋Ｐｉ，ｊ－ｋ!"Ｄｚ２
$%Ｐｉ，ｊ＝"２
Ｐ"ｚ
２
ｉ，ｊ
＝１２!!"ｚ２Ｎ
ｋ＝１&$ｋ
Ｐｉ，ｊ＋２ｋ－２Ｐｉ，ｊ＋Ｐｉ，ｊ－２ｋ!"油气地球物理２００７年７月
?６?
李胜军：波动方程的变步长有限差分数值模拟第５卷第３期?７?
若在深度ｚ０＝ｌ!ｚ的步长是双倍的，（９）式在长
度为２Ｎ的有限差分法可被使用到深度为（ｌ－Ｎ）!ｚ，（７）式可在深度ｌ!ｚ以下使用。对剩余的（Ｎ－１）个网格点，可以结合（９）和（１０）式来计算其导数。在深度（ｌ－Ｎ＋ｎ）!ｚ的网格上的导数可写作
如果Ｚ方向上的步长通过其他参数来调整，类似的方程可被推导。一般通过函数ｍ（ｋ）来产生，有
差分因子!ｋ必须根据函数ｍ（ｋ）计算，这可通过Ｆｏｒｎｂｅｒｇ（１９８８ａ）提出的一种算法来实现［９］。图４表示了方程（１２）在双倍步长区域差分算子２Ｎ＝１０的情况下有限差分参数的变化情况。显示了６个网格
列表（黑圆点位置属于精细网格循环；方框点属于双倍的步长循环点；五角星位置点的导数被计算）。在每一列用不同的差分模式来计算在不同深度的二阶导数。针对十阶精度有限差分相应的ｍ（ｋ）的值列于表１，其他阶精度的ｍ（ｋ）值可用相应的方法计算。
为保证这种差分算法的稳定性，压力必须内插在水平方向两倍步长边界上的一些点上。对参数长度２Ｎ＝１０的那些点在图５中已列出。内插计算时可使用特别精确且没有错误发生的三角内插算法。
在图２显示的网格中，可使用这种变网格算法来改变空间分辨率。在细网格区域不需要用傅立叶变换法或高阶有限差分算子，四阶或六阶精度的算法就可以满足精度的要求。
３算例分析
为验证上述算法的正确性、有效性，设计了图６
所示的低速地表地质模型。数值模拟时，震源设置在（８００ｍ，２４ｍ）处。
在内存为１．０ＧＢｙｔｅ、ＣＰＵ为奔腾ＩＶ－２．８ＧＨｚ的微型计算机上模拟试算，分别用常网格有限差分算法和变网格有限差分进行模拟运算形成单炮记录。在此基础上对４ｍ常网格步长、８ｍ常网格步长、变网格步长算法从内存需求和计算时间进行比较分析。图７（ａ）是常网格８ｍ步长模拟的单炮记录，可以看出，由于在低速地表的采样点太少，图中的频散现象比较严重，甚至影响到了倾斜界面的反射记录。
（１１）
Ｄｚ２
!"Ｐｉ，ｊ＝１!#$ｚ２Ｎ－ｎ
ｋ＝１
%!ｋＰｉ，ｊ＋ｋ－２Ｐｉ，ｊ＋Ｐｉ，ｊ－ｋ#
$＋
１２!#$ｚ２Ｎ
ｋ＝１
%!ｋＰｉ，ｊ＋２ｋ＋#$ｎ－２Ｐｉ，ｊ＋Ｐｉ，ｊ－２ｋ＋#$ｎ#$（１２）
Ｄｚ２
!"Ｐｉ，ｊ＝
１!#$ｚ２Ｎ
ｋ＝１
%!ｋＰｉ，ｊ＋ｍ（ｋ）－２Ｐｉ，ｊ＋Ｐｉ，ｊ－ｍ（ｋ）#$ｚ０
ａ
ｂ
ｃ
ｄ
ｅ
ｆ
图４过渡带的差分格式
ｍ（ｋ）ａｂｃｄｅｆｍ（０）００００００ｍ（１）１１１１１２ｍ（２）２２２２３４ｍ（３）３３３４５６ｍ（４）４４５６７８ｍ
（５）５
６
７
８
９
１０
表１图４显示的有限差分参数ｍ（ｋ）值
ｘ
ｚ
ｚ０
图５必须内插压力的网格点
１０００ｍ／ｓ
１８００ｍ／ｓ
２２００ｍ／ｓ
８００１６０００
（ｍ）
图６低速地表地质模型
０８００
１６００
（ｍ）
油气地球物理２００７年７月?８?
为压制频散，提高模拟精度，只有减小采样间隔，如用常网格４ｍ步长算法进行模拟，内存需求将会扩大１倍，以上述（１６００ｍ×１６００ｍ）的模拟区域为例，计算时间将会增大１．２倍左右。用文中所述变步长算法进行模拟，低速地表用４ｍ小步长，地表以下用８ｍ大步长，变网格步长算法的内存需求仅是常网格４ｍ步长算法的６５％左右，计算时间也缩短４５％左右。与常网格８ｍ步长算法相比，计算时间是常网格８ｍ步长的１．２倍，当然，内存需求也相应有所扩大。图７（ｂ）是变网格步长算法的模拟结果（为了对比明显，两图使用了相同的增益），可以看出，频散得到了较好的压制，分辨率也明显提高。
４结论
上述分析与理论模型试算结果表明，变网格步长差分算法是一种高效的正演模拟算法。它能较好地实现对低速地表、高速层夹层、复杂界面的数值模拟，成功地解决小网格步长差分算法内存需求量大、计算效率低等问题。在特殊地区用增大（减小）步长的模拟方法减小了内存需求，提高了计算效率；在特定位置减小步长，解决了对超薄层地质体、倾斜界面的数值模拟效果不理想的问题：解决了用常规高阶有限差分或傅立叶变换法模拟存在的问题。
该算法通过插值法来计算过渡带的波场值，在网格大小变化梯度不大、模拟精度不很高（过渡带比较窄）时，这种变网格算法是较好的。但是，如果网格调整较大、模拟精度较高，将会在过渡边界处产生一个弱的人为边界反射，同时插值边界也可能会存在稳定性问题。因此，该算法仍有一些不完善的地方。今后可以把变网格算法进一步简化，取消网格的横向变化来避免插值计算，增加算法的稳定性。
参考文献
［１］余德平，波场数值模拟技术．勘探地球物理进展，２００４，２７（１）：１６￣２１［２］Ｄ．ＫｏｓｌｏｆｆａｎｄＥ．Ｂａｙｓａｌ．ＦｏｒｗａｒｄｍｏｄｅｌｉｎｇｂｙａＦｏｕｒｉｅｒｍｅｔｈｏｄ．Ｇｅｏｐｈｙｓｉｃｓ，１９８２，４７（１０）：１４０２￣１４１２
［３］Ｍ．Ａ．Ｄａｂｌａｉｎ．Ｔｈｅａｐｐｌｉｃａｔｉｏｎｏｆｈｉｇｈ－ｏｒｄｅｒｄｉｆｆｅｒｅｎｃｉｎｇｔｏｔｈｅｓｃａｌａｒｗａｖｅｅｑｕａｔｉｏｎ．Ｇｅｏｐｈｙｓｉｃｓ，１９８６，５１（１）：５４￣６６
［４］ＣｏｒｄＪａｓｔｒａｍ，ＡｌｆｒｅｄＢｅｈｌｅ．Ａｃｏｕｓｔｉｃｍｏｄｅｌｉｎｇｏｎａｇｒｉｄｏｆｖｅｒｔｉ－ｃａｌｌｙｖａｒｙｉｎｇｓｐａｃｉｎｇ．ＧｅｏｐｈｙｓｉｃａｌＰｒｏｓｐｅｃｔｉｎｇ，１９９２，４０（２）：１５７￣１６９
［５］Ｊ．Ｅ．Ｖｉｄａｌｅ．Ｃｏｍｍｅｎｔｏｎ‘Ａｃｏｍｐａｒｉｓｏｎｏｆｆｉｎｉｔｅ－ｄｉｆｆｅｒｅｎｃｅａｎｄＦｏｕｒｉｅｒ－ｍｅｔｈｏｄｃａｌｃｕｌａｔｉｏｎｓｏｆｓｙｎｔｈｅｔｉｃｓｅｉｓｍｏｇｒａｍｓ’ｂｙＣ．Ｒ．Ｄａｕｄｔｅｔａｌ．ＢｕｌｌｅｔｉｎｏｆｔｈｅｓｅｉｓｍｏｌｏｇｉｃａｌｓｏｃｉｅｔｙｏｆＡｍｅｒｉｃａ．１９９０，８０：４９３￣４９５
［６］ＤａｎＫｏｓｌｏｆｆ，Ａ．Ｑ．Ｆｉｌｈｏ，Ｅ．Ｔｅｓｓｍｅｒ．Ｎｕｍｅｒｉｃａｌｓｏｌｕｔｉｏｎｏｆｔｈｅａｃｏｕｓｔｉｃａｎｄｅｌａｓｔｉｃｗａｖｅｅｑｕａｔｉｏｎｓｂｙａｎｅｗｒａｐｉｄｅｘｐａｎｄｉｏｎｍｅｔｈｏｄ．ＧｅｏｐｈｙｓｉｃａｌＰｒｏｓｐｅｃｔｉｎｇ．１９８９，３７（４）：３８３￣３９４
［７］孙成禹，张吉辉．完全纵波方程有限差分波场模拟．石油地球物理勘探，２００５，４０（３０）：２８９￣２９４
［８］Ｏ．Ｈｏｌｂｅｒｇ．Ｃｏｍｐｕｔａｔｉｏｎａｌａｓｐｅｃｔｓｏｆｔｈｅｃｈｏｉｃｅｏｆｏｐｅｒａｔｏｒａｎｄｓａｍｐｌｉｎｇｉｎｔｅｒｖａｌｆｏｒｎｕｍｅｒｉｃａｌｄｉｆｆｅｒｅｎｔｉａｔｉｏｎｉｎｌａｒｇｅ－ｓｃａｌｅｓｉｍ－ｕｌａｔｉｏｎｏｆｗａｖｅｐｈｅｎｏｍｅｎａ．ＧｅｏｐｈｙｓｉｃａｌＰｒｏｓｐｅｃｔｉｎｇ１９８７，３５（６）：６２５￣６５５
［９］Ｂ．Ｆｏｒｎｂｅｒｇ．Ｇｅｎｅｒａｔｉｏｎｏｆｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅｆｏｒｍｕｌａｓｏｎａｒｂｉｔｒａｒｉｌｙｓｐａｃｅｄｇｒｉｄｓ．ＭａｔｈｅｍａｔｉｃｓｏｆＣｏｍｐｕｔａｔｉｏｎ，１９８８，５１：６９９￣７０６
８００１６０００（ｍ）０．０
０．５
１．０
１．５
（ｓ）
８００１６０００（ｍ）０．０
０．５
１．０
１．５
（ｓ）
（ａ）均匀网格（ｂ）变网格
图７常网格与变网格得到的炮集记录
（下转第１３页）
彭才：动态反褶积方法研究
第５卷第３期?１３?
Ｔｉｍｅ－ＦｒｅｑｕｅｎｃｙＴｒａｎｓｆｏｒｍｓ．Ａｌｂｅｒｔａ．ＵｎｉｖｅｒｓｉｔｙｏｆＣａｌｇａｒｙ，１９９９：３～３２
［６］Ｍａｒｇｒａｖｅ，Ｇ．Ｆ，Ｌａｍｏｕｒｅｕｘ，Ｍ．Ｐ．Ｇａｂｏｒｄｅｃｏｎｖｏｌｕｔｉｏｎ：ＣＲＥＷＥＳＲｅｓｅａｒｃｈＲｅｐｏｒｔ．２００１，１３：２７８～３１０
［７］ＡｌａｎａＲ．Ｓｃｈｏｅｐｐ．ＩｍｐｒｏｖｉｎｇＳｅｉｓｍｉｃＲｅｓｏｌｕｔｉｏｎｗｉｔｈＮｏｎｓｔａｔｉｏｎａｒｙＤｅｃｏｎｖｏｌｕｔｉｏｎ．Ａｌｂｅｒｔａ．ＵｎｉｖｅｒｓｉｔｙｏｆＣａｌｇａｒｙ，１９９８：６７～８５
［８］Ｍａｒｇｒａｖｅ，Ｇ．Ｆ，Ｌａｍｏｕｒｅｕｘ，Ｍ．Ｐ．ＧａｂｏｒｄｅｃｏｎｖｏｌｕｔｉｏｎｏｆｓｅｉｓｍｉｃｄａｔａｆｏｒｓｏｕｒｃｅｗａｖｅｆｏｒｍａｎｄＱｃｏｒｒｅｃｔｉｏｎ，Ｇｅｏｐｈｙｓｉｃｓ，１９９８，７２：２１９０～２１９３
［９］Ｉｌｉｅｓｃｕ，Ｖ，Ｍａｒｇｒａｖｅ，Ｇ．Ｆ．ＧａｂｏｒｄｅｃｏｎｖｏｌｕｔｉｏｎａｐｐｌｉｅｄｔｏａＢｌａｃｋｆｏｏｔｄａｔａｓｅｔ：ＣＲＥＷＥＳＲｅｓｅａｒｃｈＲｅｐｏｒｔ．２００１，１３：２４１～２７６［１０］彭才．动态子波估计和动态反褶积方法研究［学位论文］．四川成都：西南石油大学，２００７
Ｒｅｓｅａｒｃｈｏｆｎｏｎｓｔａｔｉｏｎａｒｙｄｅｃｏｎｖｏｌｕｔｉｏｎｍｅｔｈｏｄ
ＰＥＮＧＣａｉ１，２）（１．ＧｅｏｐｈｙｓｉｃａｌＥｘｐｌｏｒａｔｉｏｎＣｏｍｐａｎｙ，ＳＰＡ，Ｃｈｅｎｇｄｏｕ６１０５０９，２．ＳｏｕｔｈｗｅｓｔＰｅｔｒｏｌｅｕｍＵｎｉｖｅｒｓｉｔｙ，６１０５００，Ｃｈｅｎｇｄｏｕ，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：Ｔｈｅｒｅａｒｅｍａｎｙｍｅｔｈｏｄｓｏｆｄｅｃｏｎｖｏｌｕｔｉｏｎａｔｐｒｅｓｅｎｔ，ｂａｓｅｄｏｎｃｏｎｖｅｎｔｉｏｎａｌｃｏｎｖｏｌｕｔｉｏｎｍｏｄｅｌ．Ｉｎｃｏｎｖｅｎｔｉｏｎａｌｃｏｎｖｏｌｕｔｉｏｎｍｏｄｅｌ，ｓｕｐｐｏｓｉｎｇｔｈｅｗａｖｅｌｅｔｎｏｔｔｏｂｅａｔｔｅｎｕａｔｅｄ，ｓｏｉｔｄｏｅｓｎ＇ｔａｃｃｏｒｄｗｉｔｈａｃｔｕａｌｉｔｙｔｈａｔｔｈｅａｍｐｌｉｔｕｄｅａｔｔｅｎｕａｔｉｏｎｇｒａｄｕａｌｌｙ．Ｔｈｅｒｅｆｏｒｅ，ｉｔ’ｓｎｅｃｅｓｓａｒｙｔｏｅｓｔａｂｌｉｓｈａｎｏｎｓｔａｔｉｏｎａｒｙｃｏｎｖｏｌｕｔｉｏｎｍｏｄｅｌｗｈｉｃｈｃａｎｄｅｓｃｒｉｂｅｔｈｅｎｏｎｓｔａｔｉｏｎａｒｙｗａｖｅｌｅｔ．Ｂａｓｅｄｏｎｔｈｅｃｏｎｖｏｌｕｔｉｏｎｍｏｄｅｌ，ｎｏｎｓｔａｔｉｏｎａｒｙｄｅｃｏｎｌｕｔｉｏｎｃａｎｉｍｐｒｏｖｅｔｈｅｒｅｓｏｌｕｔｉｏｎｏｆｓｅｉｓｍｉｃｄａｔａｗｉｔｈｏｕｔｃｏｍｐｅｎｓａｔｔｉｎｇａｍｐｌｉｔｕｄｅｏｆｓｅｉｓｍｉｃｔｒａｃｅ．Ｃｏｍｐａｒｉｎｇｔｈｅｒｅｓｕｌｔｓｏｆｓｐｉｋｉｎｇｄｅｃｏｎｖｏｌｕｔｉｏｎａｎｄｎｏｎｓｔａｔｉｏｎａｒｙｄｅｃｏｎｖｌｕｔｉｏｎｉｎｔｈｅｐａｐｅｒ，ｔｈｅａｎａｌｙｓｉｓｏｆｔｈｅｏｒｅｔｉｃａｌｍｏｄｅｌａｎｄｒｅａｌｄａｔａｃｏｍｐｕｔａｔｉｏｎｐｒｏｖｅｔｈｉｓｍｅｔｈｏｄｉｓｒｉｇｈｔ．
Ｋｅｙｗｏｒｄｓ：ｎｏｎｓｔａｔｉｏｎａｒｙｄｅｃｏｎｖｏｌｕｔｉｏｎ，ｎｏｎｓｔａｔｉｏｎａｒｙｗａｖｅｌｅｔ，ｓｅｉｓｍｉｃｄａｔａａｎｄａｎａｌｙｓｉｓｏｆｔｉｍｅ－ｆｒｅｑｕｎｃｙ
（上接第８页）
Ｗａｖｅｅｑｕａｔｉｏｎｎｕｍｅｒｉｃａｌｍｏｄｅｌｉｎｇｏｎａｇｒｉｄｏｆｖａｒｙｉｎｇｓｐａｃｉｎｇ
ＬＩＳＨｅｎｇ－ｊｕｎ（ＣｏｌｌｅｇｅｏｆＧｅｏ－ＲｅｓｏｕｒｃｅａｎｄＩｎｆｏｒｍａｔｉｏｎ，ＣｈｉｎａＵｎｉｖｅｒｓｉｔｙｏｆＰｅｔｒｏｌｅｕｍ，Ｄｏｎｇｙｉｎｇ，２５７０６１，Ｃｈｉｎａ；ＮｏｒｔｈｗｅｓｔＩｎｓｔｉｔｕｔｅｏｆＲｅｓｅａｒｃｈＩｎｓｔｉｔｕｔｅｏｆｐｅｔｒｏｌｅｕｍＥｘｐｌｏ－ｒａｔｉｏｎａｎｄＤｅｖｅｌｏｐｍｅｎｔ，ＣＮＰＣ，Ｌａｎｚｈｏｕ７３００２０，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：Ｔｈｅｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅｉｓｏｎｅｏｆｃｏｍｍｏｎｆｏｒｗａｒｄｍｏｄｅｌｉｎｇｍｅｔｈｏｄｓ，ｗｈｉｃｈｃｏｎｔａｉｎｓａｂｕｎｄａｎｔｓｅｉｓｍｉｃｉｎｆｏｒｍａ－ｔｉｏｎ，ａｎｄｉｓｅａｓｙｔｏｂｅｒｅａｌｉｚｅｄ．Ｔｒａｄｉｔｉｏｎａｌｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅｕｓｅｕｓｕａｌｌｙｕｎｉｆｏｒｍｇｒｉｄｓ，ｔｈｅｍｅｔｈｏｄｓｉｓｌａｃｋｏｆｆｌｅｘｉｂｉｌｉｔｙｉｎｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎｏｆｌｏｗ－ｓｐｅｅｄ／ｈｉｇｈ－ｓｐｅｅｄｍｅｄｉｕｍ，ｔｈｉｎｌａｙｅｒ／ｔｈｉｃｋｌａｙｅｒ．Ｉｎｔｈｉｓｐｒｅｓｅｎｔａｔｉｏｎ，ａｎａｌｇｏｒｉｔｈｍｏｆｖａｒｙｉｎｇｍｅｓｈｉｓｒｅｃｏｍｍｅｎｄｅｄ，ｗｈｉｃｈｈａｓｓｏｌｖｅｄｔｈｅｐｒｏｂ－ｌｅｍｄｅｓｃｒｉｂｅｄａｂｏｖｅ．Ｔｈｅｄｉｆｆｅｒｅｎｃｅｏｆｓｔｏｒａｇｅｒｅｑｕｉｒｅｍｅｎｔａｎｄｃｏｍｐｕｔａｔｉｏｎａｌｅｆｆｏｒｔｂｅｔｗｅｅｎｔｒａｄｉｔｉｏｎａｌｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅａｎｄｔｈｅａｌｇｏｒｉｔｈｍｏｆｖａｒｙｉｎｇｍｅｓｈｉｓｄｉｓｃｕｓｓｅｄ．Ａｎｄｔｈｅｂｏｕｎｄａｒｙｃｏｎｄｉｔｉｏｎ，ｒａｐｉｄｅｘｐａｎｓｉｏｎｍｅｔｈｏｄ（ＲＥＭ）ｆｏｒｔｈｅｔｉｍｅｉｎｔｅｇｒａｔｉｏｎｏｆｔｈｅａｌｇｏｒｉｔｈｍｏｆｖａｒｙｉｎｇｍｅｓｈａｒｅｉｎｔｒｏｄｕｃｅｄ，ｔｈｅａｄｖａｎｔａｇｅｓｏｆｔｈｉｓａｌｇｏｒｉｔｈｍａｒｅｓｕｍｍａｒｉｚｅｄ．ＫｅｙＷｏｒｄｓ：ｖａｒｙｉｎｇｓｐａｃｉｎｇ，ｂｏｕｎｄａｒｙｃｏｎｄｉｔｉｏｎ，ｃｏｍ－ｐｕｔｅｔｉｍｅ，ｒａｐｉｄｅｘｐａｎｓｉｏｎｍｅｔｈｏｄａｎｄｎｕｍｅｒｉｃａｌ
ｍｏｄｅｌｉｎｇ
叠加地震记录的相移波动方程正演模拟数值模拟实验共22页
《地震数值模拟》实验报告 一、实验题目 叠加地震记录的相移波动方程正演模拟
二、实验目的 1.掌握各向同性介质任意构造、水平层状速度结构地质模型的相移波动方程正演模拟基本理论 2.实现方法与程序编制 3.由正演记录初步分析地震信号的分辨率。 三、实验原理 1、地震波传播的波动方程 设（x,z）为空间坐标，t为时间，地震波传播速度为v(x,z)，则二位介质中任意位置、任意时刻的地震波场为p(z,x,t)：压缩波——纵波。则二维各向同性均匀介质中地震波传播的遵循声波方程为 2、傅里叶变换的微分性质 p(t)与其傅里叶变换的P(w)的关系： 3、地震波传播的相移外推公式 令速度v不随x变化，只随z变化，则利用傅里叶变换微分性质把波动方程(变换到频率-波数域，得： 4、初始条件和边界条件 按照爆炸界面理论，反射界面震源在t=0时刻同时起爆，此时刻的波场就是震源。根据不同情况，可直接使用反射系数脉冲或子波作震源。如果直接使用反射系数作震源脉冲，则初始条件可表示为： 5、边界处理
（1）边界反射问题 把实际无穷空间区域中求解波场的问题化为有穷区域求解时，左右两边使用零边界条件。物理上假设探区距Xmin与Xmax两个端点很远，在两个端点上收到的反射波很弱。但是，上述条件在实际中不能成立，造成零边界条件反而成为绝对阻止波通过的强反射面。在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。 （2）边界强反射的处理 镶边法、削波法、吸收边界都能有效消除边界强反射。 削波法就是在波场延拓过程中，没延拓一次，在其两侧均匀衰减到零，从而消除边界强反射的影响。假设横向总长度为NX，以两边Lx道吸波为例，有以下吸波公式： 四、实验内容
波动方程的变步长有限差分数值模拟
收稿日期：２００７－０３－２３；修订日期：２００７－０４－２７ 作者简介：李胜军，男，在读硕士研究生，研究方向为地震波传播理论。联系电话：（０５４６）８３９２０５５，Ｅ－ｍａｉｌ：ｈｄｐｕｌｉｓ＠１２６．ｃｏｍ，通讯地址：（２５７０６１）中国石油大学（华东）地球资信与信息学院。 ＊中国石油大学（华东）研究生创新基金资助，编号：Ｓ２００６—０６。 油气地球物理 ２００７年７月 ＰＥＴＲＯＬＥＵＭＧＥＯＰＨＹＳＩＣＳ 第５卷第３期 在地震资料采集、处理和解释中通常需要进行地震波场数值模拟：假设已知地下的地质情况，应用地震波运动学和动力学的基本原理，计算给定地质模型的地震响应。这种做法对正确认识地震波的运动学和动力学特征，以及准确分析油气藏的反射波场特征有着重要的指导意义。声波在介质中的正演模拟研究为我们精确模拟地震波在复杂介质中的传播提供了理论基础［１］。 傅立叶变换法和高阶有限差分法（ＦＤ）已成为计算声波方程空间导数的标准技术［２，３］。虽然常网格步长差分算法比较容易实现，但是它们对大部分模型都增大了不必要的计算量。例如，对存在浅层低速带的沉积盆地模型地面地震记录进行模拟时，由于低速地层阻抗小，地震波传入其中会引起较大的振幅和较长的延续时间（这与深层的高速层完全不同）。由于这些浅层低速层中地震波的波长较短、地层厚度较小，模拟时需要用小网格进行。这样，常网格步长算法就必须用小网格离散整个模型，从而增加了不必要的代价，如内存、计算量的增大。 因而，采用变网格算法将能改进有上覆低速层情况模拟结果的有效性（对地层中间有超薄夹层的情形，必须用精细网格覆盖才能精确的对地层进行模拟）。应用这种变网格算法既能实现对夹层的模拟，又能保障计算量不增加。因此这种通过函数实现在任意深度上网格步长变化的有限差分方法被 推广［４］。为了计算空间导数，在Ｘ方向用傅立叶变换法或有限差分算法，在Ｚ方向使用高阶有限差分方法。通过时间积分快速展开法（ＲＥＭ）来保障差分方法的计算精度［６］。这种差分技巧比二阶时间差分有较高的精确度且计算用时短。 １时间积分 均匀介质中的二维声波方程可用下式表示［２］ 式中：Ｐ＝Ｐ（ｘ，ｚ，ｔ），代表压力项；ｃ＝ｃ（ｘ，ｚ），代表速度；ｓ＝ｓ（ｓ，ｚ），代表震源函数；Ｌ２为差分算子。在密度!＝!（ｘ，ｚ）变化的情况下，常用的是Ｖｉｄａｌｅ给出的公式［５］ 波动方程的变步长有限差分数值模拟＊ 李胜军１，２） 孙成禹１） 张玉华１） 倪长宽１） １）中国石油大学地球资源与信息学院；２）中石油勘探开发研究院西北分院 摘要：有限差分算法是常用的正演模拟方法之一，其包含的地震信息丰富，且实现简单。传统的有限差分方法通常都采用均匀网格步长，在对含低速／高速介质、 薄层／厚层介质的模型进行波场模拟时往往缺乏稳定性。文章介绍了一种可以有效解决上述问题的变网格算法，对常规有限差分法与变网格差分算法在内存需求、计算速率等方面的差别进行了比较，对变网格差分算法中的边界条件、 时间积分的快速展开算法作了阐述，进而总结了变网格算法的优点。关键词：变步长；边界条件；计算时间；快速展开法；数值模拟 !２ Ｐ!ｔ２＝－Ｌ２Ｐ＋ｓ （１） （２） －Ｌ２ ＝ｃ ２ !２!ｘ２＋!２ !ｚ ２" # （３） （４） !２ Ｐ!ｔ ２＝－Ｌ２Ｐ!"$ －１!Ｌ２Ｐ＋ＰＬ２１!＋ｓ －Ｌ２ ＝!ｃ ２ ２ !２!ｘ２＋!２ !ｚ ２% $
声波方程数值模拟实验报告
声波方程数值模拟实验报告 一．基础理论知识 需要的已知条件包括： 1．1）震源函数 2）地层速度（波速） 3）边界条件 2．弹性波方程：?????????+??=??+??+??=??) ()()(222222 22222 222z w x w v t w t S z u x u v t u s p 声波方程的有限差分法数值模拟 对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述： )()(2222 222t S z u x u v t u +??+??=?? （4-1） (,)v x z 是介质在点（x , z ）处的纵波速度，u 为描述速度位或者压力的波场，)(t s 为震 源函数。 为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。为此，先把空间模型网格化（如图4-1所示）。 设x 、z 方向的网格间隔长度为h ?，t ?为时间采样步长，则有： h i x ?= （i 为正整数） h j z ?= (j 为正整数)t n t =? (n 为正整数) k j i u , 表示在(i,j)点，k 时刻的波场值。 将1 ,+k j i u 在(i,j)点k 时刻用Taylor 展式展开： z ?,i j 1,i j +2,i j +1,i j -2,i j -,2 i j -,2 i j +,1i j +,1 i j -1,1i j -+1,2 i j -+2,1i j -+2,2 i j -+1,2 i j ++2,2 i j ++1,1 i j +-2,1i j +-1,1i j ++2,1i j ++1,1i j --1,2i j +-2,2i j +-2, 2i j --2,1 i j --1,2i j --x ?
波动方程正演模型的研究与应用
波动方程正演模型的研究与应用 郑鸿明* 娄 兵 蒋 立 （新疆油田公司勘探开发研究院地物所） 摘要野外采集的地震数据是经过大地滤波后的畸变信号，处理的地震剖面只是间接地反映了地下构造和地质体的特征，虽然目前有很多方法和手段可以分析并提取相关的地质信息，但由于处理对波场的改造和噪声的存在以及方法本身的多解性问题降低了识别地质信息的可靠性。处理中每一步对有效信息的影响有多大，对地震属性解释的影响有多大，没有一个定量的标准，只能凭经验和认识来定性地判断。正演模型在弹性波理论指导下，遵循严格的数学公式，可以最佳模拟地下各种情况。各种处理方法和不同的处理流程所得到的结果能否符合或最佳逼近波动方程建立的数学模型，正演模型是判断处理工作合理性的良好准则。 主题词地质模型波动方程正演模型地震响应模块测试 1 引 言 随着地震勘探的不断深入，地震勘探也由构造型油气藏勘探进入精细的岩性勘探阶段，要求地震勘探能够反映地下地质体岩性变化，以及识别含油、气、水的地震响应特征，分辨薄互层、低幅度构造的能力。地球物理学家们在长期的实践中已经研究开发了很多相关的技术，虽然理论上这些方法都能够成立，这些技术应用成功的实例也很多，但也不乏有失败的教训，往往产生多解性，或与钻探的结论不符。这里除了复杂地表和复杂地下构造形成的复杂地震波场而不满足建立在简单地质模型处理理论的因素外，与处理过程对地震波场的改造也有很大关系。从地震数据的采集到最终处理的地震剖面，整个过程是一个系统工程，地下地质结构、地质体的岩性变化以及含流体的性质，对处理人员来说是看不见、摸不着的“黑匣子”，我们所看到的只是经过大地滤波后产生畸变的地震波场，如何从这个畸变的地震波场中去伪存真、恢复真实的构造形态、提取储层的相关地震属性信息，这是岩性处理的最终目标。处理中的每一步环环相扣、相互影响、相互制约，而我们对处理中的每一步产生的中间结果所应达到的标准只是凭经验、感觉进行定性判定，加入了很多人为因素，这些因素或多或少影响着我们对解释成果的正确认识。另外，处理技术发展很快，相应的地震处理软件越来越多，应用这些模块之前对各模块所起的作用以及它们所产生的结果都需要有一个定量的认识，以及验证处理流程的合理性是当前迫切需要解决的问题。究竟什么样的结果满足岩性解释的要求、什么样的结果反映的是真正地下地质体的响应、什么样的处理方法满足保振幅处理和地震属性分析的应用等等一系列问题，这都是当前岩性处理中迫切需要解决的主要问题。它直接关联着处理成果的真伪及后续解释的可靠性，关联着勘探的投资风险。 随着计算机运算能力发展迅猛，特别是微机群的出现，为波动方程算法提供了硬件环境，开展此项技术的研究与应用已成为可能。此次模型的设计全面考虑了地表和地下的典型地质特征并将这些特征容入到模型中，真实模拟了实际地质结构。应用该地质模型正演叠前炮集的地震响应。 2 模型的建立 模型分物理模型和数学模型两种，目前的物理模型只能做非常简单的模拟，只有用数学模型才能模拟各种复杂的地质现象。20世纪70年代，美国哥伦比亚大学在郭宗汾
2007射线追踪与波动方程正演模拟方法对比研究
４７ 科技资讯 科技资讯 ＳＣＩＥＮＣＥ ＆ ＴＥＣＨＮＯＬＯＧＹ ＩＮＦＯＲＭＡＴＩＯＮ２００７ ＮＯ．１２ ＳＣＩＥＮＣＥ ＆ ＴＥＣＨＮＯＬＯＧＹ ＩＮＦＯＲＭＡＴＩＯＮ 工 业 技 术 地震正演模拟作为反演解释的反过程，是验证解释成果的有效手段，进行必要可靠的正演模拟可以有效的监控反演解释。地震学一般可以分为几何地震学和物理地震学，在几何地震学中进行的正演模拟方法就是我们通常所说的射线追踪法，射线追踪法是在合成记录时用地震子波和界面或地质体的反射系数进行反褶积运算，即。运算的最大特点是说明了地震波传播的运动学特征。而在物理地震学中应用波动方程法合成的地震记录是通过求解波动方程的数值解来模拟地震波场的。在波动方程合成的地震记录中不单保持了地震波传播运动学特征，还说明了地震波传播的动力学特征。本文将分别用射线追踪和波动方程的方法合成地震记录。 １ 基于射线追踪的合成地震响应 射线追踪法的主要理论基础是，在高频近 射线追踪与波动方程正演模拟方法对比研究 王志美 畅永刚 （长江大学油气资源与勘探技术教育部重点实验室 湖北荆州 ４３４０２３） 摘 要：地震学一般可以分为几何地震学和物理地震学，几何地震学中进行正演模拟方法就是射线追踪法，射追踪法是在合成记录时用地震子波和界面或地质体的反射系数进行反褶积运算，即。运算的最大特点就是说明了地震波传播的运动学特征。而在物理地震学中的波动方程法合成的地震记录是通过求解波动方程的数值解来模拟地震波场。在波动方程合成的地震记录中不单保持了地震波传播 运动学特征外，还说明了地震波传播的动力学特征。本文将分别用射线追踪和波动方程的方法合成地震记录。关键词：射线追踪 波动方程 正演模拟 中图分类号：Ｐ３１５文献标识码： Ａ 文章编号：１６７２－３７９１（２００７）０４（ｃ）－００４７－００ 图１ 射线追踪正演模拟（１） 图 ２ 逐段迭长示意图 图 ３ 射线追踪正演模拟（２） 图４ 波动方程正演模拟结果 似条件下，地震波的主能量沿射线轨迹传播。基于这种认识，运用惠更斯原理和费马原理来重建射线路径，并利用程函方程来计算射线的旅行时。在旅行时计算中应用有限差分等方法，以获得快速的解。射线法的主要优点是概念明确，显示直观，运算方便，适应性强；其缺陷是应用有一定限制条件，计算结果在一定程度上是近似的，对于复杂构造进行两点三维射线追踪往往比较麻烦。为了计算波沿射线的旅行时和波的传播路径，叙述如下。 如图１所示，首先给出连接Ｓ（激发点）和Ｒ（接收点）之间的初始射线路径射线的振幅变化，首先必须知道地震波在实际地层中传播的射线路径。 由于地震波在整条路径上满足同一个射线参数，因此射线路径上任意连续三点也将满足同一个参数，而三点间的射线表现形式为Ｓｎｅｌｌ定律。按照Ｓｎｅｌｌ 定律，可导出一个求 取中间点的一阶近似公式。当前后两点位于界面两边时，中间点为透射点，所求路径为透射路径；当前后两点位于界面的同一边时，中间点为反射点，所求路径为反射路径。为此，可以从任一端点出发，连续地选取三点，通过一阶近似公式进行逐段迭代取中间点，利用新求出的点代替原来的点，然后以一点的跨跃作为步长，顺序地逐段迭代下去，直到另一端点。这样，新计算出的中间点和两个端点就构成了一次迭代射线路径，如图２中所示。如果整条射线路径上校正量的范数之和满足一定的精度要求，则认为射线追踪过程结束，否则从追踪出的射线路径开始，继续重复上述过程，直到满足精度要求为止。最后一次追踪到的中间点和两个端点，构成整条射线路径。图３基于多层倾斜界面模型通过射线追踪正演模拟地震响应。从模拟结果可以直观的看出基于几何地震学的原理正演模拟结果只能反映地震波的几何传播路径。在实际的工程设计中通过正演模拟可以在地表确定地下观测范围，节约设备提高工程效率，但不能反映 物理地震学中的地震属性，例如振幅，频率和相位等。更不能反映地震波的动力学特征。 ２ 波动方程的合成地震响应 ２．１ 波动方程的建立 非均匀介质的声波方程： （１） （２） 可由对连续介质方程（１）式的两端对时间求导，并利用欧拉方程推得： （３） 其中：Ｐ是波数，Ｖ是质点振动的速度向量，ρ是密度，ｃ是波速，ρ和ｃ是随着空间参数χ和ｚ变化的，这里ρ给定为常数，只有ｃ 是地质模型的控制参数。χ和Ｚ分别是在地面水平距离和深度。这样（３）式就可以变为： （４） 其中：ｃ＝ν （χ，ｚ）；（４）式即是所求的弹性波动方程。 ２．２ 数值计算及稳定性 求解弹性波动方程的方法有多种，付立叶变换法是对弹性波动方程的波场进行付立叶变换，优点是运算速度快。克希霍夫积分法是基于均匀模型，利用格林函数公式计算曲面积分，求出空间波场值，但这种方法不能适应
地震波数值模拟方法研究综述.
地震波数值模拟方法研究综述 在地学领域，对于许多地球物理问题，人们已经得到了它应遵循的基本方程(常微分方程或偏微分方程)和相应的定解条件，但能用解析方法求得精确解的只是少数方程性质比较简单，且几何形状相当规则的问题。对于大多数问题，由于方程的非线性性质，或由于求解区域的几何形状比较复杂，则不能得到解析解。这类问题的解决通常有两种途径。一是引入简化假设，将方程和几何边界简化为能够处理的情况，从而得到问题在简化状态下的解答。但这种方法只是在有限的情况下是可行的，过多的简化可能导致很大的误差甚至错误的解答。因此人们多年来寻找和发展了另一种求解方法——数值模拟方法。 地震数值模拟(SeismicNumericalModeling)是地震勘探和地震学的基础，同时也是地震反演的基础。所谓地震数值模拟，就是在假定地下介质结构模型和相应的物理参数已知的情况下，模拟研究地震波在地下各种介质中的传播规律，并计算在地面或地下各观测点所观测到的数值地震记录的一种地震模拟方法。地震波场数值模拟是研究复杂地区地震资料采集、处理和解释的有效辅助手段，这种地震数值模拟方法已经在地震勘探和天然地震领域中得到广泛应用。 地震数值模拟的发展非常迅速，现在已经有各种各样的地震数值模拟方法在地震勘探和地震学中得到广泛而有效
的应用。这些地震波场数值模拟方法可以归纳为三大类，即几何射线法、积分方程法和波动方程法。波动方程数值模拟方法实质上是求解地震波动方程，因此模拟的地震波场包含了地震波传播的所有信息，但其计算速度相对于几何射线法要慢。几何射线法也就是射线追踪法，属于几何地震学方法，由于它将地震波波动理论简化为射线理论，主要考虑的是地震波传播的运动学特征，缺少地震波的动力学信息，因此该方法计算速度快。因为波动方程模拟包含了丰富的波动信息，为研究地震波的传播机理和复杂地层的解释提供了更多的佐证，所以波动方程数值模拟方法一直在地震模拟中占有重要地位。 1地震波数值模拟的理论基础 地震波数值模拟是在已知地下介质结构的情况下，研究地震波在地下各种介质中传播规律的一种地震模拟方法，其理论基础就是表征地震波在地下各种介质中传播的地震波传播理论。上述三类地震波数值模拟方法相应的地震波传播理论的数学物理表达方式不尽相同。射线追踪法是建立在以射线理论为基础的波动方程高频近似理论基础上的，其数学表形式为程函方程和传输方程。积分方程法是建立在以惠更斯原理为基础的波叠加原理基础上的，其数学表达形式为波动方程的格林函数域积分方程表达式和边界积分方程表达式。波
基于GPU的波动方程正演模拟的实现
基于GPU的波动方程正演模拟的实现 袁崇鑫;邓飞 【期刊名称】《电脑知识与技术》 【年(卷),期】2014(000)018 【摘要】随着计算机技术的发展，使得波动方程正演由理论研究应用到实际地震勘探中成为了可能。而有限差分技术作为地震波场模拟的一种有效数值方法，它具有实现简单，速度快，从而被广泛应用正演计算密集的波形正反演中。地震波正演的计算量大，通过CPU来计算地震波正演模拟严重影响整体运算效率，GPU通用计算技术的产生及其在内的数据并行性有望改变这一状况。该文主要研究波动方程正演在GPU上的模拟实现。%With the development of computer technology, the wave equation forward by the application of theory to real seismic exploration as possible. The finite-difference seismic wave field simulation technology as an effective numerical methods, it has a simple, fast, and thus is widely used computationally intensive forward modeling and inversion of the waveform. Computationally intensive seismic forward modeling of seismic waves through the CPU to calculate the forward modeling seriously affect the over-all operational efficiency, GPU general computing technologies, including the generation and data parallelism is expected to change this situation. This paper studies the wave equation forward simulation on the GPU. 【总页数】6页(4333-4337,4340)
声波波动方程正演模拟程序总结
声波波动方程正演模拟程序 程序介绍： 第一部分：加载震源，此处选用雷克子波当作震源。 编写震源程序后，我将输出的数据复制，然后我用excel做成了图片，以检验程序编写是否正确。以下为雷克子波公式部分的程序： for(it=0;it 模型构建与试算： 1、我首先建立了一个均匀介质模型，首先利用不同时间，进行了数值模拟，得到波场快照如图所示： 100ms 200ms 300ms 此处，纵波速度为v=3000m/s。模型大小为200×200，空间采样间隔为dx=dz=10m。采用30Hz的雷克子波作为震源子波，时间采样间隔为1ms，图中可以看出，波场快照中的同相轴是圆形的，说明在均匀各向同性介质中，点源激发的波前面是一个圆，这与理论也是吻合的。并且随着时间的增大，波前面的面积逐渐增大，说明地震波从震源中心向外传播。 2、我在建立的均匀模型的基础上，改变差分算子的精度，分别采用2阶、6阶、12阶精度进行试算。时间统一采用300ms的时候。得到的波长快照如下： 2阶精度6阶精度12阶精度 第44卷 第2期 煤田地质与勘探 Vol. 44 No.2 2016年4月 COAL GEOLOGY & EXPLORA TION Apr . 2016 收稿日期: 2014-11-06 作者简介: 智敏(1984—)，男，江苏盐城人，助理研究员，从事地震资料处理方法研究. E-mail: zhimin@https://www.360docs.net/doc/463355290.html, 引用格式: 智敏. 声波方程隐式高阶有限差分数值模拟[J]. 煤田地质与勘探，2016，44(2)：106–111. ZHI Min. High-order implicit finite difference numerical simulation of acoustic wave equation[J]. Coal Geology & Exploration ，2016，44(2)：106–111. 文章编号: 1001-1986(2016)02-0106-06 声波方程隐式高阶有限差分数值模拟 智 敏 (中煤科工集团西安研究院有限公司，陕西 西安 710077) 摘要: 推导了声波方程空间二阶导数的隐式求解公式及差分系数的求解方法，讨论了该方法的数值频散特征。利用该方法分别对均匀介质及Marmousi 模型进行了数值模拟，将其结果与传统的显式差分格式的模拟结果进行了对比分析。结果表明：该方法较传统的显式求解方法具有更低的数值频散、更高的计算精度。 关 键 词：声波方程；数值模拟；隐式差分；数值频散 中图分类号：P631 文献标识码：A DOI: 10.3969/j.issn.1001-1986.2016.02.019 High-order implicit finite difference numerical simulation of acoustic wave equation ZHI Min (Xi'an Research Institute , China Coal Technology and Engineering Group Corp , Xi'an , 710077, China ) Abstract: The paper derived implicit solution method for solving differential equations and coefficient of acoustic space of the second derivative of the equation, discussed numerical frequency dispersion characteristics of this method. Using this method, Marmousi model and homogeneous medium were simulated numerically. The results were analyzed and compared with the result obtained by traditional explicit difference scheme. The result showed that this method was lower in numerical dispersion and more accurate in calculation. Key words: acoustic wave equation; numerical simulation; implicit difference; numerical frequency dispersion 声波方程数值模拟已广泛应用于地震正演、逆时偏移、全波形反演等各个方面，如何进一步提高数值模拟精度、减小频散和计算量一直是有限差分算法研究的核心问题[1-2]。目前，相关文献给出的空间二阶导数计算大多采用显式差分格式，这种方法通常具有较强的频散现象[3-4]。而为了有效地降低数值频散，通常采用减小网格尺寸、提高差分阶数的方法，这种办法虽然提高了数值模拟的精度，但也大幅增加了计算量和内存空间[5]。Boris 和Book 将求解流体力学方程的通量校正传输方法(FCT)应用到波动方程数值模拟中，该方法通过差分运算降低高波数端能量压制了频散噪声，但不利于提高模拟结果的分辨率[3-7]。 采用隐式差分格式求解波场的空间二阶导数较显式法具有更高的精度。前人研究表明：隐式差分格式的(2N +2)阶差分算子的精度与显式差分算子(4N +2)阶的计算精度相当[8] 。因此，笔者研究了空间二阶导数的隐式求解方法，理论分析及模拟结果均表明：该方法较显式求解方法具有更小的空间频散及更高的模拟精度。 1 方法原理 1.1 基本方程 二维声波方程如式(1)所示： 2222 2 22 1p p p x y v t ???+=??? (1) 式中 p 为地震波场；v 为纵波速度。如果只考虑x 方向，根据波场传播理论可有：()()i e k x p x x p x ?+?=，因此根据欧拉公式可推出下式： ()()()() () () 2224sin /2p x p x x p x x p x k x p x x x x δδ+?+-?--?= = ?? (2) 式中 22() p x x δδ为差分算子；k 为波数；x ?为空间采样间隔。 为了提高差分精度，可以取： 224 2224 x p p p q b x x x x δδδδ?=≈-?? (3) 式中 b 为一待定常数。由于 1 1,11x x x ≈-<<+，所以 地震波波动方程数值模拟方法 地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶变换法、有限元法和有限差分法等。 克希霍夫积分法引入射线追踪过程，本质上是波动方程积分解的一个数值计算，在某种程度上相当于绕射叠加。该方法计算速度较快，但由于射线追踪中存在着诸如焦散、多重路径等问题，故其一般只能适合于较简单的模型，难以模拟复杂地层的波场信息。 傅里叶变换法是利用空间的全部信息对波场函数进行三角函数插值，能更加精确地模拟地震波的传播规律，同时，利用快速傅里叶变换(FFT)进行计算，还可以提高运算效率，其主要优点是精度高，占用内存小，但缺点是计算速度较慢，对模型的适用性差，尤其是不适应于速度横向变化剧烈的模型. 波动方程有限元法的做法是：将变分法用于单元分析，得到单元矩阵，然后将单元矩阵总体求和得到总体矩阵，最后求解总体矩阵得到波动方程的数值解；其主要优点是理论上可适宜于任意地质体形态的模型，保证复杂地层形态模拟的逼真性，达到很高的计算精度，但有限元法的主要问题是占用内存和运算量均较大，不适用于大规模模拟，因此该方法在地震波勘探中尚未得到广泛地应用。。 相对于上述几种方法，有限差分法是一种更为快速有效的方法。虽然其精度比不上有限元法，但因其具有计算速度快，占用内存较小的优点，在地震学界受到广泛的重视与应用。 声波方程的有限差分法数值模拟 对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述： )()(2222222t S z u x u v t u +??+??=?? （4-1） (,)v x z 是介质在点（x , z ）处的纵波速度，u 为描述速度位或者压力的波场，)(t s 为震源函数。 为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。为此，先把空间模型网格化（如图4-1所示）。 设x 、z 方向的网格间隔长度为h ?，t ?为时间采样步长，则有： z ?,i j 1,i j +2,i j +1,i j -声波方程隐式高阶有限差分数值模拟
地震波波动方程数值模拟方法