声波方程正演模拟1
声波方程数值模拟实验报告
声波方程数值模拟实验报告.基础理论知识需要的已知条件包括：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。
声波方程有限差分正演
题目：使用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 ∆≤即为所求得的稳定性条件。
声波波动方程正演模拟分析研究
Δt
c
在数值计算中,生 成 的 强 边 界 反 射 会 对 中 心 波
其中:
σ0 =l
og
1
R
3
x
δ
2
(
11)
3Vp
,
R 是 理 论 反 射 系 数;
δ
2δ
是 PML 的厚 度;Vp 为 速 度;在 此 基 础 上 可 推 演 由
PML 边界条件进行交错网格的有限差分格式。
理论上,
PML 法对各种入射角和频率下的地震
弥散,产生数值频散现象 [9],严重影响正演模拟的精
度。为了减轻数值频散,通常可采用以下方式:① 调
整恰当的时间和空 间 离 散 步 长,尤 其 空 间 步 长 不 宜
过大过小,而时间步长相对越小越好,但会受到实际
计算效率的限制。 ② 通 过 提 高 差 分 阶 数,其 中 提 高
空间差分阶数易实 现 且 能 有 效 降 低 频 散,但 随 着 空
2023 年 7 月
第 13 期 总第 527 期
Ju
l
y2023
No.
13 To
t
a
lNo.
527
内 蒙 古 科 技 与 经 济
I
nne
r Mongo
l
i
aSc
i
enc
eTe
chno
l
ogy & Ec
onomy
声波波动方程正演模拟分析研究
朱晓洁
(中国石化胜利油田分公司海洋采油厂,山东 东营 257237)
波的吸收效果良好,吸 收 效 率 强 于 传 统 的 吸 收 边 界
条件法 [14]。因此在计算区域加入吸 收 边 界 后,可 以
正演计算声波时差的理论曲线
谢谢大家！
LT
2l
O2
O1
O
深度记录点：
O O1 O2 2
可消除深度误差
声系2
优点：减小井眼扩径的影响 消除深度误差 减小测量误差 分辨率比双发双收要高
缺点：操作不便
声系3：单发四收
优点：分辨率高 可得到四种不同分辨率的时差
缺点：深度误差 扩径影响
参数设定：
砂岩声速：4500m/s 泥岩声速：1800m/s 泥浆声速：1500m/s 井径：0.25m
声系1
不可识别 0.5m厚及 以下地层， 地层厚度大 于0.5m时 半幅点对应 地层厚度
声系2
可识别 0.5m厚的 地层，地层 厚度大于 0.5m时半 幅点对应地 层厚度
声系3
可识别0.15m 厚的地层，地 层厚度大于 0.15m时，半 幅点对应地层 厚度，扩径处 的上下界面时 差会偏大和偏 小
正演计算声波时差理论曲线
声系1：双发双收
t增大 T1
R1
扩 径
T1和T2交替发射声脉冲，
井 段
分别测量时差 t1和t 2。
t1 最终记录的声波时差为：
R2
t 2
未 扩 径
t t1 t2 2
井
段
T2
优点
T1
（基本）消除了扩径的影响
可消除（减小）深度误差R1Biblioteka o 盲 实际传播路径中点：
R2
o o
区
o o o 2
缺点
T2
分辨率降低
对低速地层会出现“盲区”
R1 R2
UT LT
B A
声系2：
位置A：测LT与R1，R2之间
波动方程正演模拟边界条件的比较分析
波动方程正演模拟边界条件的比较分析付小波;韩超;原健龙;余嘉顺【摘要】通过数值模拟研究了透明边界、Clayton-Engquist边界和完全匹配层边界的吸收效果,得出如下结论:在反射角和频率相同的情况下,完全匹配层边界条件效果最好,Clayton-Engquist边界效果次之,而透明边界条件的效果最差.以边界条件对100 Hz模型边界垂直反射的吸收效果来衡量,完全匹配层边界条件与Clayton-Engquist边界条件的效果分别是透明边界条件的16.5倍和3.5倍.在＜40 Hz的低频范围内,或者在反射角＞65°的情况下,Clayton-Engquist边界相对透明边界的吸收效果相对优势显著变弱.而完全匹配层边界的吸收效果则在150 Hz频率范围内和75°反射角范围内始终保持稳定的相对优势.【期刊名称】《成都理工大学学报（自然科学版）》【年(卷),期】2015(042)004【总页数】8页(P492-499)【关键词】波动方程;正演模拟;有限差分;边界条件【作者】付小波;韩超;原健龙;余嘉顺【作者单位】成都理工大学地球物理学院,成都610059;成都理工大学地球物理学院,成都610059;成都理工大学地球物理学院,成都610059;成都理工大学地球物理学院,成都610059;新西兰皇家地质与核科学研究所,惠灵顿【正文语种】中文【中图分类】P631.4边界条件是地震波数值模拟方法技术中的一项重要内容。
许多专家学者从不同角度提出多种构造边界条件的方法。
1977年，Clayton与Engquist［1］根据旁轴近似理论（Claerbout［2］；Claerbout与Johnson［3］），提出利用一系列不同近似精度的单程波动方程来吸收模型边界的反射能量，称作Clayton-Engquist （下文采用简略记号CE来表达）边界条件。
1978年，Reynolds通过对波动方程的分解得到了透明边界条件［4］（下文采用简略记号TBC来表达）。
声波波动方程正演模拟程序总结
声波波动方程正演模拟程序程序介绍：第一部分：加载震源，此处选用雷克子波当作震源。
编写震源程序后，我将输出的数据复制，然后我用excel做成了图片，以检验程序编写是否正确。
以下为雷克子波公式部分的程序：for(it=0;it<Nt;it++){t1=it*dt;t2=t1-t0;source[it]=(1.0-2.0*pow(PI*fm*t1,2.0))*exp(-pow(PI*fm*t1,2.0));fprintf(fp,"%.8f %.8f\n",t1,source[it]);}此处，为了成图完整，我用的是t2，而不是t1，也就是把雷克子波向右移动了一段距离，使主要部分都显示出来。
（频率采用的是30hz)从图中可以看出程序是正确的，符合理论上雷克子波的波形。
第二部分：主程序，编写声波正演模拟算子。
首先定义了各种变量，然后指定震源位置，选择权系数，给速度赋值，然后是差分算子的编写，这是主要部分，最后再进行时间转换，即把n-1时刻的值给n时刻，把n时刻的值给n+1时刻。
此处，我编写的是均匀介质声波方程规则网格的正演模拟程序，时间导数采用二阶中心差分、空间导数为2N阶差分精度，网格大小为200*200，总时间为400。
第三部分：这一部分就是记录文件。
首先记录Un文件，然后记录record文件。
模型构建与试算：1、我首先建立了一个均匀介质模型，首先利用不同时间，进行了数值模拟，得到波场快照如图所示：100ms 200ms 300ms此处，纵波速度为v=3000m/s。
模型大小为200×200，空间采样间隔为dx=dz=10m。
采用30Hz的雷克子波作为震源子波，时间采样间隔为1ms，图中可以看出，波场快照中的同相轴是圆形的，说明在均匀各向同性介质中，点源激发的波前面是一个圆，这与理论也是吻合的。
并且随着时间的增大，波前面的面积逐渐增大，说明地震波从震源中心向外传播。
声音的波动方程
声音的波动方程声音是一种能够通过空气、水等介质传导的物理现象，而声音的波动方程就是描述声波在介质中传播时的数学公式。
以下将从几个步骤来阐述声音的波动方程。
第一步：介质的振动声波是由介质分子的振动引起的，当声波在介质中传播时，它们会引起介质的周期性振动。
因此，声波可以被视为机械波，与大多数其他类型的波一样，它是由波的振幅、频率和波长三个要素确定的。
第二步：波的传播速度声波的传播速度取决于介质的密度、弹性模量和介质的压缩性等因素。
根据拉普拉斯原理，和波源和接收器之间的距离有关，声波的传播速度可以写成一个公式：v=fλ其中v是声波的传播速度，f是声波的频率，λ是声波的波长。
第三步：声波的压强变化声波的传播是通过介质压强变化的方式来实现的。
当声波通过介质时，它们会引起介质的压缩和膨胀。
这导致压强在空气中产生变化，使空气分子在颤动。
通过这种方式，声波在空气中传播。
第四步：声波的波动方程声波的波动方程可以用下列偏微分方程表示：∇²p(x,y,z)-1/v²*(∂²p(x,y,z)/∂t²)=0其中∇²是拉普拉斯算子，p（x、y、z，t）是压强（即声波的幅度）的空间和时间变化，v是声波的传播速度。
在坐标系中，x、y、z表示空间位置变量，t表示时间变量。
因此，这个方程可以解释为“空间中压强的二阶时间倒数等于时间中压强的拉普拉斯算子除以速度的平方”。
结论声音的波动方程是根据物理原理得出的，在声波传播的所有过程中都起到了关键作用。
通常情况下，声波在介质中的传播速度、波长、频率和波幅等特性由声音的波动方程计算得出。
因此，声音的波动方程是研究声波性质和声学的重要基础。
声波方程正演模拟共46页
t
x y z
vx 1 (P )
t x
vy 1 (P )
t y
vz 1 (P )
t z
11
二、波动方程类型及其局限性
能够描述且只能描述纵波的传播规律，包括 直达波、反射波、透射波、折射波等，但不能描 述转换波传播规律。
需要的已知条件包括： 1）震源函数 2）地层速度/密度 3）边界条件
对于二维速度-深度模型，地下介质中地震波 的传播规律可以近似地用声波方程描述：
2u t 2
v
2
(
2u x2
2u z 2
)
S (t )
(4-1)
v(x, z) 是介质在点（x , z）处的纵波速度，
u 为描述速度位或者压力的波场，
s(t) 为震源函数。
23
空间模型网格化（如图4-1所示）：
i 2, j2
12
2、弹性波方程：
2u t 2
v
p
2
(
2u x 2
2u z 2 ) S (t)
2
w
t 2
vs
2
(
2w x 2
2w z 2 )
vx
xx
xz
t
x
z
vz xz zz
t
x
z
xx ( 2u) vx vz
t
x
z
zz vx ( 2u) vz
t
x
z
xz u vz u vx
内容提纲ห้องสมุดไป่ตู้
一、地震勘探基本原理 二、波动方程类型及其 局限性 三、数值算法类型及其局限性 四、声波方程的有限差分法数值模拟
1
一、地震勘探基本原理
＊
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
傅里叶变换法是利用空间的全部信息对波场函数 进行三角函数插值，能更加精确地模拟地震波的传播 规律，同时，利用快速傅里叶变换(FFT)进行计算，还 可以提高运算效率，其主要优点是精度高，占用内存 小，但缺点是计算速度较慢，对模型的适用性差，尤 其是不适应于速度横向变化剧烈的模型。
20
波动方程有限元法的做法是：将变分法用于单元分 析，得到单元矩阵，然后将单元矩阵总体求和得到总体 矩阵，最后求解总体矩阵得到波动方程的数值解；其主 要优点是理论上可适宜于任意地质体形态的模型，保证 复杂地层形态模拟的逼真性，达到很高的计算精度，但 有限元法的主要问题是占用内存和运算量均较大，不适 用于大规模模拟，因此该方法在地震波勘探中尚未得到 广泛地应用。
13
2、弹性波方程：
2u t 2
v
p
2
(
2u x 2
2u z 2 ) S (t)
2
w
t 2
vs
2
(
2w x 2
2w z 2 )
vx
xx
xz
t
x
ห้องสมุดไป่ตู้
z
vz xz zz
t
x
z
xx ( 2u) vx vz
声波方程数值模拟
――地球物理场论 基础Ⅰ期末作业（1）
任课教师： 宋鹏
1
内容提纲
一、地震勘探基本原理 二、波动方程类型及其 局限性 三、数值算法类型及其局限性 四、声波方程的有限差分法数值模拟
2
一、地震勘探基本原理
＊
▽ ▽▽▽ ▽
x
t
3
同相轴为双曲线，即反射波的时距曲线为双 曲线，反射波一个同相轴可带来一个地层的信息。 实际地下介质非常复杂，所得到的炮集记录也包含 更多的地下信息。实际的炮集记录见图1-1和1-2。
10
内容提纲
一、地震勘探基本原理 二、波动方程类型及其 局限性 三、数值算法类型及其局限性 四、声波方程的有限差分法数值模拟
11
二、波动方程类型及其局限性
1、声波方程：
二阶标量声波方程： 一阶压力-速度方程组：
2 p t 2
v2
(
2 p x2
2 p z 2
)
S(t)
P C 2 (vx vy vz )
对于二维速度-深度模型，地下介质中地震波 的传播规律可以近似地用声波方程描述：
2u t 2
v
2
(
2u x2
2u z 2
)
S (t)
(4-1)
v(x, z) 是介质在点（x , z）处的纵波速度，
u 为描述速度位或者压力的波场，
s(t) 为震源函数。
24
空间模型网格化（如图4-1所示）：
4
图1-1 陆上某区实际地震记录
5
图1-2 海上某区实际地震记录
6
广义的地震反演即是从地震炮集记录出发，经 过复杂的去噪、速度分析以及偏移成像处理等手段 得到反映地下的地质结构的地震剖面。实际的地震 剖面见图1-3和1-4。
7
图1-3 陆上某区地震剖面
8
图1-4 海上某区地震剖面
9
地震波场模拟即地震正演，是指已知模型结构， 通过物理或数值计算的方法模拟该地质结构下的地 震波的传播，最终合成地震记录，也可以认为其是 野外数据采集过程的室内再现。物理模拟花费昂贵， 人们一般采用比较经济的数值模拟技术。地震波场 数值模拟是在给定数学模型（如弹性波方程，声波 方程等）、震源和地下几何界面、物性参数（岩层 密度、速度等）情况下，研究弹性波或声波的传播 规律。
t
x
z
zz vx ( 2u) vz
t
x
z
xz u vz u vx
t
x
z 14
能够描述纵、横波的传播规律，包括直达波、反 射波、透射波、折射波以及转换波等。
需要的已知条件包括： 1）震源函数 2）地层速度或根据方程的类型需要提 供的地层的其它弹性参数 3）边界条件
17
内容提纲
一、地震勘探基本原理 二、波动方程类型及其 局限性 三、数值算法类型及其局限性 四、声波方程的有限差分法数值模拟
18
三、数值算法类型及其优缺点
地震波波动方程数值模拟方法主要包括克希霍夫积 分法、傅里叶变换法、有限元法和有限差分法等。
克希霍夫积分法引入射线追踪过程，本质上是波动 方程积分解的一个数值计算，在某种程度上相当于绕射 叠加。该方法计算速度较快，但由于射线追踪中存在着 诸如焦散、多重路径等问题，故其一般只能适合于较简 单的模型，难以模拟复杂地层的波场信息。
t
x y z
vx 1 (P )
t x
vy 1 (P )
t y
vz 1 (P )
t z
12
二、波动方程类型及其局限性
能够描述且只能描述纵波的传播规律，包括 直达波、反射波、透射波、折射波等，但不能描 述转换波传播规律。
需要的已知条件包括： 1）震源函数 2）地层速度/密度 3）边界条件
21
相对于上述几种方法，有限差分法是一种更为快速 有效的方法。虽然其精度比不上有限元法，但因其具 有计算速度快，占用内存较小的优点，在地震学界受 到广泛的重视与应用。
22
内容提纲
一、地震勘探基本原理 二、波动方程类型及其 局限性 三、数值算法类型及其局限性 四、声波方程的有限差分法数值模拟
23
四、声波方程的有限差分法数值模拟
15
3、粘声波/弹性波方程 前面讨论的是理想弹性介质，波在其中传播时，
没有能量的损耗，介质中应力和应变关系严格遵循 胡克定律（这种理想介质称虎克固体），但波在实 际介质中传播时，是有能量损耗的，这就是所谓的 弹性波吸收。波在传播过程中，实际介质的不同部 位之间会出现某种摩擦力，称为内摩擦力或粘滞力。 这种力导致机械能向其他形式能量转换，最终转化 为热能消耗掉。
16
在地震勘探中，地震波传播的实际介质是十分 复杂的。在一定条件下，即震源作用时间短，作用 力微小，地球介质可以看作完全弹性模型，但随着 地震勘探技术的发展，勘探精度要求提高，面临复 杂地质目标时，要求地震勘探采用更加符合实际的 介质模型进行研究。粘弹性介质模型更符合实际。
但是到目前为止，在地震资料反演处理中应用 最多的还是声波方程，弹性波以及粘弹性波方程的 应用还只是停留在模拟层次上。