声波方程有限差分正演
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
题目：使用Ricker 子波，刚性边界条件，并且初值为零，在均匀各向同性介质条件下，利用交错网格法求解一阶二维声波方程数值解。
解：
一阶二维声波方程：
22222221z
P
x P t P c ∂∂+∂∂=∂∂ (1)
将其分解为：
21P c t P
x P z x z x z V V x z V t
V t ∂∂∂⎧=+⎪∂∂∂⎪
∂∂⎪=⎨
∂∂⎪∂∂⎪=⎪∂∂⎩
(2)
对分解后的声波方程进行离散，可得到：
112211,-1,,,1
22[]N n n n n m i m j i m j xi j xi j m t V
V
c P P h +-+---=∆=+-∑ 112
2
11,1,,,12
2
[]N
n n n n m i j m i j m zi j zi j m t V
V c P P h +-++--
-=∆=+-∑ 11112
12222,,m 1,,,,11
[]N
n n n n n n i j
i j
m
xi j xi m j zi j m zi j m m tc P
P c
V
V
V
V
h
+
+
+
+
+++-+--=∆=+
-+-∑
h z x =∆=∆
针对公式(1)，使用二阶中心差商公式：
2
P(,,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 x
c 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 x
c 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 变换：
222
2exp()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 -∆-+∆⎡⎤
+⎢⎥∆-+-∆∆⎢⎥-∆-+∆⎢⎥
⎢⎥⎣⎦
2
22
222
sin sin 22sin (2
x z k x k z
w t
t c h
∆∆+∆=∆） (5) 公式(5)右端必须满足下列条件：
2
2222
sin sin 220(x z k x k z
t c h
∆∆+≤∆≤）1 取x k 和z k 最大值，即=x x k π∆，z =k z π∆，则有：2222
0t c h
≤∆≤1
因此tc ∆≤
即为所求得的稳定性条件。 在程序试算中，选中相关参数如下：
dt 0.001()10()3000/20010030N m s x z m v m s x z m
sx sz m f Hz
=⎧⎪
∆=∆=⎪⎪=⎪
==⎪⎨
⎪⎪==⎪
=⎪⎪
⎩时间网格空间网格（速度）（模型范围）t=1s(采样长度)（震源位置）（主频）=15(空间精度)
因为v 3000*0.0013t ∆==<
=，满足稳定性条件。 波长v 3000
=
100m 30f λ==，空间采样间隔h=10m ，一个波长内的空间采样点数 100
m=
1010
h
λ
=
=，基本满足网格色散条件。 程序运行输出P(,,)x z t 波场快照：
(a) 200ms (b) 300ms
(c) 400ms (d) 500ms
(e) 600ms (f) 700ms
(g) 800ms (h) 900ms
图1 波场快照
主要程序代码附录：
震源Ricker子波函数
float f(float t1, float f0)
{
float t00=1/f0,y;
y=((1.0-2.0*pow(pi*f0*(t1-t00),2))*exp(-pow(pi*f0*(t1-t00),2)));
return y;
}
计算交错网格有限差分系数
void Cal_1D_FdCoff(float *c1,int n)
{
float **A = new float * [n];