声波方程有限差分法数值模拟
声波方程数值模拟实验报告
声波方程数值模拟实验报告.基础理论知识需要的已知条件包括：1）震源函数地层速度（波速） 边界条件.2.2.2苇=v 2(諾二2)S(t)一 t :x :zv （x, Z ）是介质在点（x , z ）处的纵波速度，u 为描述速度位或者压力的波场，s （t ）为震 源函数。
为求式（4-1）的数值解，必须将此式 离散化，即用有限差分来逼近导数，用差商代替 微商。
为此，先把空间模型网格化（如图4-1所示）。
设x 、z 方向的网格间隔长度为h ，厶t 为时间采样步长，则有:X 二i^h（i 为正整数） z = j.）h（j 为正整数）^n ：t（n 为正整数）u ：j 表示在（i,j ）点，k 时刻的波场值。
k1. 2)3)>2u 2f cu〒=V p ^~2 t :x22w 2 w 厂二 V （—T :t ；x 声波方程的有限差分法数值模拟 对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述:2. 弹性波方程:二b s (t) 一 z22_w )；z 2)(4-1)将u i,j在（i,j）点k时刻用Taylor展式展开：k 1 kjUu i,j 7,j ■—ct将U i k j 」在(i,j)点k 时刻用Taylor 展式展开:k k4 k k5 k[U i2jui 2,j] 3[ui 4,j ui ・1,j] —?U i,j }S(t)*、(i -i °)**「(j - j °)(4-7)式中v(i, j)为介质速度的空间离散值，：h 是空间离散步长，=t 为时间离散步长，s(k)为震源函数，关于 s(k) 一般使用一个理论的雷克型子波代替，即:上式中，t 为时间，f 为中心频率，一般取为20-40HZ ， 为控制频带宽度的参数,(4-2)k 1k ；'UUi,jUi,2*「7 2 ；:t 2(4-3)将上两式相加，略去高阶小量, 整理得(i,j)点k 时刻的二阶时间微商为:2k 1kk J；：u u i,j -2u i,j u i,j.:t 242(4-4)对于空间微分，采用四阶精度差分格式,(以X 方向为例)即将U i*；j 、*Tj 分别在(i,j)点k 时刻展开到四阶小量，消除四阶小量并解出二阶微分得:-2~ u 11 k k 4 k k 5 k{—T7[u i 二 j +u id2,j ]+：[u i4j +U i*,j ] —=u i,j } 12 3 2:x L X 2(4-5)同理可得:2；=u 1.1 kk4 kk5 k一 2 人 2{一 石[Ui,j ，+Ui,j~2]+：[Ui,j 」+Ui,j^]—；U i,j}:z-z12 32(4-6)这就实现了用网个点波场值的差商代替了偏微分方程的微商，将上三个式子代入 (4-1)式中得:k 1 Ui,j kk 4= 2u i,j —Ui,jV :-1 1 k k 4 k k—{p [u=j U i/ yij —j：h 25 k H-U i,j }h 2 {12s (t )二 e(-2 f / )2t2cos2 二 ft(4-8)般取3-5。
有限差分法地震波传播数值模拟
2k
2 x
ω2
⎟⎟⎠⎞2
⎪⎫ ⎬ ⎪⎭
∂3P ∂z∂t 2
−
v2 4
∂3P ∂x2∂z
+
3v 4
∂3P ∂x2∂t
−
1 v
∂3P ∂t 3
=
0
Elastic Wave: (Bottom Boundary)
Utt = α 2U xx + β 2U zz + (α 2 − β 2 )Wxz Wtt = β 2Wxx + α 2Wzz + (α 2 − β 2 )U xz
2
MM 4
⎤ ⎥ ⎥
⎢⎢⎡CC12((MM
) )
⎤ ⎥ ⎥
M = M
6
⎥ ⎥ ⎥
⎢⎢C3(M
)
⎥ ⎥
⎢M⎥
2M ⎥
M⎦
⎢⎣CM(M ) ⎥⎦
⎡1⎤ ⎢⎢0⎥⎥ ⎢0⎥ ⎢⎢ M ⎥⎥ ⎢⎣0⎥⎦
⎡1
⎢ ⎢
13
⎢ 15
⎢ ⎢
M
⎢⎣12N −1
3 33 35
M 32 N −1
5 53 55
M 52 N −1
1992，1994年，Tessmer et al.在模拟二维以及三维不规则 地表面波时，同样也是使用的如上吸收边界。
采取的措施
z 高阶差分解法--提高计算精度，减小数值频散 z 采用基于特征分析方法得到的吸收边界条件
数值频散问题------高阶差分解法
声波：
∑ ∂2 f
∂x 2
=1 Δx 2
M
Cm(M )[ f (x + mΔx) − 2 f (x) +
m=1
f (x − mΔx)] + o(Δx2M )
声波方程有限差分正演
题目：使用Ricker 子波，刚性边界条件，并且初值为零，在均匀各向同性介质条件下，利用交错网格法求解一阶二维声波方程数值解。
解：一阶二维声波方程：22222221zPx P t P c ∂∂+∂∂=∂∂ (1)将其分解为：21P c t Px P z x z x z V V x z V tV t ∂∂∂⎧=+⎪∂∂∂⎪∂∂⎪=⎨∂∂⎪∂∂⎪=⎪∂∂⎩(2)对分解后的声波方程进行离散，可得到：112211,-1,,,122[]N n n n n m i m j i m j xi j xi j m t VVc P P h +-+---=∆=+-∑ 112211,1,,,122[]Nn n n n m i j m i j m zi j zi j m t VV c P P h +-++---=∆=+-∑ 1111212222,,m 1,,,,11[]Nn n n n n n i ji jmxi j xi m j zi j m zi j m m tc PP cVVVVh+++++++-+--=∆=+-+-∑h z x =∆=∆针对公式(1)，使用二阶中心差商公式：2P(,,1)2(,,)(,,1)i j n P i j n P i j n t +-+-∆222(1,,)2(,,)(1,,)(,1,)2(,,)(,1,)P i j n P i j n P i j n xc P i j n P i j n P i j n z +-+-⎧⎫+⎪⎪⎪⎪∆=⎨⎬+-+-⎪⎪⎪⎪⎩∆⎭(3)变形：P(,,1)=2(,,)(,,1)i j n P i j n P i j n +--2222(1,,)2(,,)(1,,)t (,1,)2(,,)(,1,)P i j n P i j n P i j n xc P i j n P i j n P i j n z +-+-⎧⎫+⎪⎪⎪⎪∆+∆⎨⎬+-+-⎪⎪⎪⎪⎩∆⎭(4)对离散格式作时间和空间三重Fourier 变换：0P(,,)(,,)x z i j n P k k w ↔ ，0P(,,1)(,,)*exp()x z i j n P k k w iw t +↔∆0P(1,,)(,,)*exp(k )x z x i j n P k k w i x +↔-∆，0z P(,1,)(,,)*exp(k )x z i j n P k k w i z +↔-∆对公式(4)进行Fourier 变换：2222exp()2exp()h exp()2()exp()2exp()h x x z z ik x ik x iw t iw t t c ik z ik z -∆-+∆⎡⎤+⎢⎥∆=--∆+∆⎢⎥-∆-+∆⎢⎥⎢⎥⎣⎦2222exp()2exp()h exp()2()=exp()2exp()h x x z z ik x ik x iw t iw t t c ik z ik z -∆-+∆⎡⎤+⎢⎥∆-+-∆∆⎢⎥-∆-+∆⎢⎥⎢⎥⎣⎦222222sin sin 22sin (2x z k x k zw tt c h∆∆+∆=∆） (5) 公式(5)右端必须满足下列条件：22222sin sin 220(x z k x k zt c h∆∆+≤∆≤）1 取x k 和z k 最大值，即=x x k π∆，z =k z π∆，则有：22220t c h≤∆≤1因此tc ∆≤即为所求得的稳定性条件。
二维声波方程交错网格有限差分数值模拟研究
在对介质模型进行离散化处理的过程中，网格
是一种常用手段。对波动方程进行网格离散，可以
利用交错网格的差分形式。交错网格就是把速度和
应力分配到两套不同的网格中，这样可以使速度、应
力得到很好的耦合［２］。利用交错网格有限差分法对
一阶速度—应力波动方程进行求解时，应力、速度等
分量在模型交错网格节点中的位置分布如图 １所示。
－Ｕ ］ ｋ＋１／２ ｉ，ｊ－ｎ＋１
其中 ｘ＝ｉΔｘ，ｚ＝ｊΔｚ，ｔ＝ｋΔｔ，ｉ、ｊ、ｋ分别表示空 间和时间网格点。Ｕｋｉ，＋ｊ１／２，Ｗｋｉ＋＋１１／／２２，ｊ＋１／２，Ｐｋｉ，＋ｊ＋１１／２／２，Ｑｋｉ，ｊ＋１／２ 和 Ｓｋｉ＋１／２，ｊ分别是速度 Ｖｘ、Ｖｚ与应力 σｘｘ、σｚｚ、σｘｚ的离 散值。
数的问题转化为求解网格节点上的差分方程组的问
题，得到数值解。在波动方程网格离散化的过程中，
可以利用交错网格的差分形式。
１．１ 均匀各向同性介质二维声波方程
均匀各向同性介质二维声波方程可表示为：
２ｕ（ｘ，ｘｙ２，ｚ，ｔ）＋２ｕ（ｘ，ｚｙ２，ｚ，ｔ）＝
ｖ２（１ｘ，ｚ）２ｕ（ｘ，ｔｙ２，ｚ，ｔ）
（１）
应力 Ｐｘｘ： Ｐｋｉ＋＋１１／２，ｊ＝Ｐｋｉ＋１／２，ｊ＋Ｃ１１ΔΔｘｔｎΣＮ＝１ＣＮｎ［Ｕｋｉ＋＋ｎ１，／ｊ２－Ｕｉｋ－＋ｎ１＋／２１，ｊ］＋
Ｃ１３ΔΔｚｔｎΣＮ＝１ＣＮｎ［Ｗｉｋ＋＋１１／／２２，ｊ＋ｎ－１／２
－Ｗ ］ ｋ＋１／２ ｉ＋１／２，ｊ－ｎ１／２
应力 Ｑｚｚ： Ｑｉｋ＋＋１１／２，ｊ＝Ｑｋｉ＋１／２，ｊ＋Ｃ１３ΔΔｘｔｎΣＮ＝１Ｃ（ｎＮ）［Ｕｉｋ＋＋ｎ１，／ｊ２－Ｕｉｋ－＋ｎ１＋／２１，ｊ］＋
Ｃ３３ΔΔｚｔｎΣＮ＝１Ｃ（ｎＮ）［Ｗｋｉ＋＋１１／／２２，ｊ＋ｎ－１／２－Ｗｉｋ＋＋１１／／２２，ｊ－ｎ＋１／２］
声波方程有限差分数值模拟的变网格步长算法
声波方程有限差分数值模拟的变网格步长算法声波方程是描述声波在介质中传播的方程，在进行数值模拟时需要使用有限差分方法来近似求解。
有限差分方法将连续的空间和时间离散化，通过在离散的网格节点上计算声压场的数值来模拟声波的传播过程。
在有限差分数值模拟中，选择合适的网格步长对数值结果的精度和计算效率具有重要影响。
变网格步长算法是一种能够根据需要在不同区域自适应地调整网格步长的技术。
在声波方程有限差分数值模拟中，声波会在不同介质中以不同的速度传播，因此在网格步长选择上需要考虑介质的变化。
当声波在介质变化剧烈的区域传播时，使用较小的网格步长可以提高模拟的精度。
而在介质变化缓和的区域，使用较大的网格步长能够减少计算量。
变网格步长算法的基本思想是通过对声波传播区域进行划分，然后根据介质的变化情况调整不同区域的网格步长。
具体步骤如下：1.对传播区域进行划分：根据介质的变化情况，将声波传播区域划分为多个区域，每个区域具有不同的网格步长。
2.确定初始网格步长：在初始时，可以根据经验或者初步模拟结果选择一个合适的网格步长。
3.计算声波传播：在每个区域内分别进行声波方程的有限差分数值模拟，使用当前区域的网格步长进行计算。
4.判断误差与精度：通过计算得到的声压场数值和预期精度进行比较，如果达到要求，则结束模拟；否则，进行下一步。
5.调整网格步长：根据当前区域的模拟误差情况，调整该区域的网格步长，使其更适应声波在当前区域的传播特性。
可以根据误差大小和梯度等因素进行调整。
6.重复模拟：根据调整后的网格步长，重新进行声波方程的有限差分数值模拟，并返回步骤4通过以上步骤，可以实现声波方程有限差分数值模拟的变网格步长算法。
这种算法能够根据介质的变化情况自适应地调整网格步长，提高模拟的精度和计算效率。
然而，变网格步长算法的实现并不简单，需要合理设置划分区域和调整步长的策略，并在计算过程中不断优化和调整，才能达到较好的数值模拟效果。
声波有限差分数值模拟
声波方程的离散化
时间离散化
将时间轴划分为一系列离散的时 间点，以便在每个时间点上求解 声波方程。
空间离散化
将连续的空间区域划分为一系列 离散的网格点，以便在每个网格 点上求解声波方程。
声波方程的有限差分近似
声波速度的近似
使用有限差分法近似声波的速度，以 便在离散化的空间和时间网格上求解 声波方程。
调试技巧
使用断点、单步执行和变量监视等技巧，定位和 修复代码中的错误。
04 声波有限差分法的应用
在地震波模拟中的应用
总结词
地震波模拟是声波有限差分法的重要应用领域之一，通过模拟地震波在地壳中的传播过程，可以用于地震勘探、 地震工程和地震安全性评估等领域。
详细描述
声波有限差分法能够模拟复杂的地震波场，包括反射、折射、散射等现象，为地震勘探提供精确的数据。在地震 工程中，该方法可用于评估结构的抗震性能和设计地震防护措施。此外，在地震安全性评估中，声波有限差分法 可以模拟地震对地下设施的影响，为设施的安全设计和运营提供依据。
谱方法
将谱方法与有限差分模拟 相结合，提高数值计算的 稳定性和精度。
无网格方法
探索与无网格方法（如移 动最小二乘法、径向基函 数等）的结合，实现更加 灵活的离散方式。
THANKS FOR WATC限差分法需要将连续的空间离散化为网格，这可能导致 结果的精度和可靠性受到网格大小和形状的影响。
边界条件处理
在处理具有复杂边界条件的声波问题时，有限差分法可能 需要采用特殊的边界条件处理方法，增加了模拟的复杂性 和难度。
数值色散误差
有限差分法在模拟声波传播过程中会产生数值色散误差， 即模拟结果与真实解之间的误差随模拟步数的增加而增大， 这可能影响模拟结果的精度和可靠性。
超声波在各向同性固体中传播的数值模拟_肖开丰
收稿日期:2007-08-09 修回日期:2007-08-17 第25卷 第9期计 算 机 仿 真2008年9月 文章编号:1006-9348(2008)09-0292-03超声波在各向同性固体中传播的数值模拟肖开丰,宋文爱(中北大学信息与通信工程学院,山西太原030051)摘要:介绍以应用动量守恒和应力—应变关系为各向同性弹性介质进行微变形建立超声波传播方程,并利用一种有限元方法和差分法,分析超声波(p 波)在各向同性介质中传播,及在裂缝上的散射结果。
并且通过计算机程序模拟二维超声波的传播,散射的数值模拟与几何理论结果是相符合的。
因此超声检测仿真软件是可以预测超声检测过程中的波形。
超声探头的建模发展趋势是利用F D M 和F E M 方法研究超声检测的仿真软件,分析结果将广泛应用到检测过程中。
关键词:有限元;超声波;波传播中图分类号:T P 391.9 文献标识码:BN u m e r i c a l S i m u l a t i o no f t h e P r o p a g a t i o n o f U l t r a s o n i cWa v e i na nI s o t r o p i c S o l i dX I A OK a i -f e n g ,S O N GWe n -a i(S c h o o l o f I n f o r m a t i o n a n d C o m m u n i c a t i o n E n g i n e e r i n g ,N o r t h U n i v e r s i t y o f C h i n a ,T a i y u a n S h a n x i 030051,C h i n a )A B S T R A C T :T h i s p a p e r p r e s e n t s t h e a p p l i c a t i o n o f t h e e q u a t i o n s o f c o n s e r v a t i o n o f m o m e n t u ma n d t h e s t r e s s -s t r a i n r e l a t i o n s f o r a n i s o t r o p i c e l a s t i c m e d i u mu n d e r g o i n g i n f i n i t e s i m a l d e f o r m a t i o nt o e s t a b l i s ha n u l t r a s o n i c w a v e p r o p a g a -t i o ne q u a t i o n ,a n du s e st h ef i n i t ee l e m e n t m e t h o d(F E M)a n dt h ed i f f e r e n c em e t h o df o r a n a l y s i i n gu l t r a s o n i c (P-w a v e )w a v e p r o p a g a t i o ni ni s o t r o p i c m e d i a ,a n dt h ec r a c k s s c a t t e r i n gr e s u l t s .A n dt h et w o -d i m e n s i o n a l u l t r a s o n i c w a v e p r o p a g a t i o n h a s b e e ns i m u l a t e db y c o m p u t e r p r o g r a m a n d t h e s c a t t e r e d w a v e i s w e l l e x p l a i n e db ya g e o m e t r i c a l t h e o r y .T h e r e f o r e u l t r a s o n i c t e s t i n g s i m u l a t i o n s o f t w a r e c a n p r e d i c t t h e c o u r s e o f u l t r a s o n i c t e s t i n g .T h e m o d e l s o f u l -t r a s o n i c p r o b e s w e r e d e v e l o p e d ,a n dt h e y w e r e c o m b i n e dw i t h t h e F E M a n d t h e F D M t o d e v e l o p a c o m p l e t e s i m u l a t o r o f u l t r a s o n i c t e s t i n g (U T ).T h e s e f e a t u r e s e n a b l e t h e m o d e l i n g o f a w i d e r a n g e o f u l t r a s o n i c t e s t i n g .K E Y WO R D S :F i n i t e e l e m e n t ;U l t r a s o n i c w a v e ;P r o p a g a t i o n1 引言固体中超声的散射是声学的基本问题之一,同时是超生无损检测的核心问题。
声波方程有限差分数值模拟的变网格步长算法
声波方程有限差分数值模拟的变网格步长算法声波方程有限差分数值模拟是一种常用的声波传播模拟方法，可以在计算机上通过数值计算求解声波传播的过程。
在进行这种数值模拟时，常常需要选择合适的网格步长，以保证计算结果的准确性和计算效率。
本文将介绍一种变网格步长算法，用于优化声波方程有限差分数值模拟的计算。
声波方程可以用下面的形式表示：∂^2p/∂t^2=c^2∇^2p其中p是声场变量，t是时间，c是声速，∇^2是Laplace算子。
为了将声波方程用有限差分方法进行离散化计算，我们需要将空间和时间分别离散化。
首先，将空间离散化为网格，在每个网格点上计算声场的值。
其次，将时间离散化为离散的时间步长，通过迭代计算不同时间步长上的声场分布。
为了保证计算结果的准确性，网格步长应当满足Nyquist采样定理的要求。
即网格步长应小于声波的最小波长的一半。
根据声波方程的性质，我们可以通过声速和最高频率来估计声波的最小波长。
然后，我们可以根据最小波长来选择合适的网格步长。
然而，在实际的声波传播计算中，声场的变化往往不是均匀的。
有些区域的声场变化较大，而其他区域的声场变化较小。
如果我们在整个计算区域都采用较小的网格步长，将会造成计算资源的浪费。
因此，需要一种方法能够根据声场的变化情况来自适应地调整网格步长。
变网格步长算法就是一种能够根据声场变化情况自动调整网格步长的算法。
其基本思想是根据声场在不同网格上的变化率来决定每个网格上的网格步长。
具体的算法步骤如下：1.初始化：选择一个合适的初始网格步长。
通常可以选择根据声波的最小波长来确定。
2.计算网格步长：在每个时间步长上，对于每个网格点，计算其周围网格点上的声场变化率。
常用的方法是计算声场在三个相邻时间步长上的差分值，然后取绝对值并求平均。
根据声场变化率，调整当前网格点上的网格步长。
变化率大的网格点应该有更小的网格步长，而变化率小的网格点则可以有更大的网格步长。
3.更新声场：根据调整后的网格步长，更新所有网格点上的声场值。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
}
main() {
FILE *fp; int i,j,k; float v[I][J],Dt[I][J];
//定义文件
double **u1,**u2,**u3,**u4,**u5,w[KMAX];
通过以上的计算，可以将 k=130 和 k=240 时波场的数据进行导出，用 Matlab 绘出图像。其中，在 Matlab 里用到的成像程序如下（这里只列举第一张波前快照的程序，第二张的与之相同）： clear all; load wavefront1.dat;
m=wavefront1(1); n=wavefront1(2); len=length(wavefront1)-2; nn=len/m; if n~=nn; disp('The data file length error!'); return; end
u k 1 i, j
2uik, j
uik,
1 j
v2t 2 h 2
{
1 12
[
uk i2,
j
uk i2,
j
]
4 3
[uik1,
j
uk i 1,
j
]
5 2
uik,
j
}
v 2t 2 h 2
{
1 12
[
uik,
j2
uik, j2 ]
4 3
[uik,
j
1
主要集中在前部，大致上 t>35 的时候 S(t)都为零值， 故在绘制 S(t)图像时 k 取 100，来比较清晰得呈现雷 克子波。本实验 不涉及吸 收边界条件 的使用。这样
声波方程数值模 拟所需的（ 1）震源函 数（2）地层
速度（波速）（3）边界条件均定义完全。
为求（1）式的数值解，必须将此式离散化（包 括时间离散和空间离散，这里 dt 取 0.002，dh 取 4， 即用有限差分来逼近导数，用差商代替微商。为此，
图4
图5
通过调节增益（x(j)=(i-1)*dx+wavefront((i-1)*n+j+2)*12*dx，dx 前面的 12），对波场
值放大一定倍数，可以使能量小的增大，使成像清晰。从下图中可清楚看到波前的图像，数层
的地震波随时间的推移由震源开始向外扩散。k=130 时，大约是时间为 260ms 时的图像，此时
关键词：声波方程 数值模拟 地震波 有限差分 引言：地震波场模拟即地震正演，是指已知模型结构，通过物理或数值计算的
方法模拟该地质结构下的地震波的传播，最终合成地震记录，也可以认为其是野 外数据采集过程的室内再现。物理模拟花费昂贵，人们一般采用比较经济的数值 模拟技术。地震波场数值模拟是在给定数学模型（如弹性波方程，声波方程等）、 震源和地下几何界面、物性参数（岩层密度、速度等）情况下，研究弹性波或声 波的传播规律。地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶 变换法、有限元法和有限差分法等。相对于上述几种方法，有限差分法是一种更 为快速有效的方法。虽然其精度比不上有限元法，但因其具有计算速度快，占用 内存较小的优点，在地震学界受到广泛的重视与应用。本实验给定大小两个模型， 对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模 型中间，分别记录两个时刻的波前快照（即该时刻区域内所有网格点的波场值）。
由上面的程序可得震源离散后的图像，如下图：
图2
小模型波前快照
选用一个小模型，即 200*200 的范围，视为完全弹性模型，为一层介质，在介质中取速 度为定值，本实验中取为 2500，将震源放在模型中央，分别记录两个时刻的波前快照，（即 区域内所有网格点的波场值）。第一时刻为地震波还未传播到边界上的某时刻，第二时刻为 地震波已经传播到边界上的某时刻，体会其人工边界反射。
别对应哪些波。
实验内容：
对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：
2u t 2
v
2
(2u x 2Fra bibliotek2u z 2
)
S (t )
（1）
其中 v(x,z)是介质在点(x , z)处的纵波速度，本实验中 v 为常数 2500，u 为描述速度位或者压力的波场，
s(t) 为震源函数，本实验选用雷克子波（零相位子波）
e s(t) (2f / )2t2 cos 2ft
上式中，t 为时间，fm 为中心频率，一般取为 20-40HZ，本实验取 20，γ为控制频带宽度的参数， 一般取 3-5,这里取 4，在实际计算过程中，需把此震 源函数离散，参 与波场计 算。由图像 可知子波能量
图3
根据网格设计可知，网格间距为 10m，则 200*200 的矩阵实际上大致是个 2000*2000 平方米的区域， 这样震源在（100，100）即中心点，即在地质剖面上，宽度为 1000m，深度为 1000m，假设介质均匀各向 同性，且边界为自由界面，所以在介质中速度 v 以一常数 2500m/s 的传播速度传播，应该在 400ms 左右时 到达边界。
uik, j1 ]
5 2
uik,
j
}
s(t) * (i i0 ) * ( j j0 )
（2）
一般来讲，差分时差分 格式阶数越高， 得到的波场值精 度越高，但稳定性 会受影响，差分 格式的稳定 性不仅与差分格式本身有 关，也随着差分 精度的不同而不同 ，而且与网格步 长之比的大小有 关，值得指出
u1[i][j]=0;u2[i][j]=0;u3[i][j]=0; //对速度，u1，u2，u3 赋初值
}
for(k=0;k<KMAX;k++)
{
w[k]=exp(-pow(2*PI*FM/R,2)*pow(k*DT,2))*cos(2*PI*FM*k*DT);
}
//波场值计算
for(i=0;i<I;i++)
������������ 为 Nyquist 频率，一般取为主频的两倍，G 为每个波长所占的网格点数，时间、空间为两阶差分的 情况 G 取 8，而时间、空间为四阶差分的情况 G 取 4。
则由上面两式得 v∈[1000，3061.8]
由上面定义的数值，用 Matlab 将震源函数 s(t) 成图，下面是成图的源程序：
的地震波的振幅愈来愈小。
以下为小模型波前快照的源程序： /********************************** *声波方程正演模拟——小模型波前快照* ***********************************/ #include <stdio.h> #include <math.h> #include <stdlib.h> #define PI 3.14159265 #define FM 20 #define R 4 #define KMAX 1000 #define DT 0.002 #define I 200 #define J 200 #define SI 100 #define SJ 100 #define DH 10
第一时刻为地震波还未传播到边界上的某时刻，第二时刻为地震波已经传播到边
界上的某时刻，体会其人工边界反射；对于大模型，定义为水平层状速度模型（至
少两层）；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的
波前快照，体会地震波在两种介质的分界面上传播规律；二是合成一个地震记录，
即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分
m[i]=(double*)malloc(y*sizeof(double));
return m;
}
void exchange(double **u1,double **u2,double **u3) //循环交换赋值 {
int i,j; for(i=0;i<I;i++)
for(j=0;j<J;j++) u1[i][j]=u2[i][j];
3
vmax * t / h
的是差分格式的稳定性与微分方程无关，本实验所采用的（2）式为条件稳定，稳定条件：
8，
（2）式由差分方程近似替代微分方程所得，所以如果空间和时间采样间隔不当，就会产生频散现象，导致
波形畸变，甚至派生出多个同相轴，为减少频散，Δh 需满足 Dablain 的经验公式： h vmin /(GfN )
声波方程有限差分法数值模拟
学 院：海洋地球科学学院 专业：地球信息科学与技术 年级：2012 级 姓 名：王昊 指 导教师：宋鹏
声波方程有限差分法数值模拟
2012 级 地球信息科学与技术 12040032037 王昊
摘要：本实验应用声波方程作为正演模拟的波动方程，将所提供震源函数离散
后绘图，并将给定两个二维速度-深度模型（一个小模型；一个大模型），绘出图 形。通过模拟地震波在介质中的传播，理解实际勘探中地震波在地层中的传播规 律，在模拟水平层状速度模型中，体会地震波在两种介质分界面的传播规律，并 能够从地震记录中识别出反射波，透射波，多次波，折射波和绕射波。并通过模 拟人工合成的地震记录，体会地震勘探基本原理和方法，验证地震波传播能量波 形变化趋势。
在小模型中，可将介质看做是单一介质，即在整个介质中，波的传播速度是相同的，在 这里设地震波速度为 2500m/s，以下是 Matlab 成图的源程序：
clf for x=1:200
f(x)=2500; end plot(f);axis([1,200,1,3000]);title('小模型的速度-深度模型'); figure; x=[0 0 200 200]; y=[0 200 200 0]; patch(x,y,'g'); title('小模型'); axis ij; 用 matlab 编译后得图像，见下图：