有限差分法波动方程正演模型制作方法研究.ppt
傅里叶有限差分法三维波动方程正演模拟
ｓｇ ｉｃ ｔ ｅ ｕ ｅ ｏ ｕ ａｉｎ ｏ ｔｃ ｕ ｅ ｙ ｅ ｒｒｃ ｒｅ ｔ ｎ．Ｃｏ ａｅ ｔ ｈ ｎｔ － ｉ ｅｅ ｃ ｔｏ ｉｎｆ ａ ｌ ｒｄ ｃ ｓｃ ｍｐ ｔｔ ａ ｃ ｓ ａ ｓｄ ｂ ｒｏ ｏ ｃｉ ｉｎ ｙ ｏｌ ｏ ｍｐ ｄ ｗ ｈ ｔ ｅ ｆ ｉ ｄｆ ｒｎ ｅ ｍｅ ｄ ｒ ｉ ｉ ｅ ｈ ｂ ｈ ｚ ｒ－ ｆｓｔ ｒ ｃ ｒ ｓ ｎ ｄ ｈ ｔ ｐ ｏ ｌｓ ａ ｅ ｏ － ｙ ｔｅ ｅ ｏ ｏ ｅ ｅ ｏｄ ａ ｓ ｏ ｒｆ ｅ ｂ ｄ ｎ ３ Ｄ ｍｏ ｉｅ Ｆ ｅ ｃ ｄ ｅ ．ｈ ｒ ｐ ｓ ｄ ｉ ｓ ｄｆ ｄ ｒｎ ｈ ｍｏ １ ｔｅ ｐ ｏ ｏ ｅ ｍｅ ｏ ａ ｉ ｈ ｔｄ ｃ ｎ ｐ ｐ ｒ ｄ ｌｔｅ ｐ ｍａ ｒｆ ｃｅ ｖ ｓａ ｄ ｉ ｃ ｒ ｔｒｃｉｅ ｉ ｏ ｏ ｏ ｒ ｅｌ ｍｏ ｅ ｈ ｒ ￣ ｅ ｅ ｔｄ ｗａ ｅ ｙ ｉ ｌ ｎ ｓ ｍｕ ｈ ｍｏ ｅ ａｔａ ｔｖ ｎ ｂ ｔ ｃ ｍｐｕａｉｎ ｏｔａ ｄ ｓｏ ａ ｅ ｈ ｔｔ ａ ｃ ｓ ｏｌ ｎ ｔｒｇ
Ｃ ｉｅｅ ， ０ ７ ０ ６ ：８ ４— １６ ｈｎ ｓ ） ２ ０ ，５ （ ） １ ５ ８２
傅 里 叶有 限差 分 法 三 维 波 动 方 程 正 演模 拟
张金 海 ， 卫 民 ， 连 锋 ， 振 兴 王 赵 姚
中 国科 学 院地 质 与地 球 物 理 研 究 所 ， 京 １０ ２ 北 ００９
Ｍ ｏ ｌｎ － ｓ ａ ａ ｖ ｓ ｕ ｉ ｇ ｔ ｕ ｉ ｒ ｆ ｉｅ ｄｉ ｅ ｅ ｃ ｅ ｈ ｄｅｉ ｇ ３ Ｄ ｃ ｌ ｒ ｗａ ｅ ｓｎ ｈｅ Ｆｏ ｒ ｅ ＇ ｔ － ｆ ｒ ｎ ｅ ｍ ｔ ｏ ｍ ｄ
Ｚ ＡＮ ｉ－ ｉ Ｗ Ａ ｉ ｎ， Ｈ ｉｎ Ｆ ｎ ， ＡＯ Ｚ ｅ － ｉｇ Ｈ Ｇ Ｊ Ｈａ ， ＮＧ Ｗｅ－ ｎ Ｍｉ Ｚ ＡＯ Ｌ ａ － ｅ ｇ Ｙ ｈ ｎ Ｘ ｎ
傅里叶有限差分法三维波动方程正演模拟
( x , y , z + Δz ;ω) . P′
( 3)
最后 ,将瞬时波场 P ( x , y , z + Δz ;ω) 代入 ( 4 ) 式进 行有限差分校正 , 9 P ( x , y , z + Δz ;ω) 9z
1856
2 2
地 球 物 理 学 报 ( Chinese J . Geophys. )
[22 ] [19 ] [13 ,14 ]
可以进一步改进相应的算法 . 基于 ADIPI 格式的 FFD 法既消除了分裂误差 , 又不影响 x 方向和 y 方向的精度 , 而且保持了双向 分裂 FFD 法的高效性 , 因此它可以胜任高陡倾角 、 强横向变化介质的三维正演模拟 ,成为精度高 、 效率 高、 适应性强的三维单程波正演方法 .
2 基于 ADIPI 格式的 FFD 方法
211 三维 FFD 算子
横向非均匀介质中单程波方程的三维 FFD 算 子为
[16 ,17 ]
9 9 2 + 2 9 x 9 y 0 Δs + kz ≈ k z + ω 2 2 9 9 1 + a 2 + 2 9x 9y
b
2
2
,
( 1)
其中 ω 是角频率 , k z =
Zhang J H ,Wang W M , Zhao L F ,et al. Modeling 3 2D scalar waves using the Fourier finite2difference method. Chinese J . Geophys . ( in Chinese) , 2007 , 50 (6) :1854～1862
1 引 言
波动方程差分方法初步(PPT文档)
U
n1 j
(1
2 2 )U
n j
2
U
n j 1
U
n j 1
U
n1 j
U
n 0
(n
)
n,
U
n J
(n )
n,
n0
U
0 j
fj,
U
1 j
fj
gj
2
2
f j1 2 f j f j1 ,
0 jJ
h
U
1 j
2(1
2
)U
0 j
2
U
0 j 1
U
0 j 1
U
1 j
初始速度的离散
一、简单处理
初始位移
U
0 j
f
(
jh)
f
0 j
初始速度 u(xj ,tk ) u(x j ,tk ) u(x j ,tk ) O( )
t
u(x j ,tk ) u(x j ,tk ) u(x j ,tk )
U
0 j 1
U
1 j
2U
1 j
2(1
2
)U
0 j
2
U
0 j 1
U
0 j 1
2 g j
2f j 2 g j 2 f j1 2 f j f j1
U
1 j
fj
gj
第五章 有限差分法 知识讲解课件
的 m=4，即此表对应差商的精度是四阶的。从这些表可以看出，一般地说，随着
差分阶数的增大和对应差商精度的提高，差分表达式所包含的项数将增多。
表 5-1
j
n0 1 2 34
1 -1
aj 1
2 1 -2 1
3 -1 3 -3 1
4 1 -4 6 -4 1
表 5-3 j
n0 1 2345 aj
1 -3 4 -1 2 2 -5 4 -1 3 -5 18 -24 14 -3 4 3 -14 26 -24 11 -2
依此类推，任何阶差分都可由其低一阶的差分再作一阶差分得到。例如 n 阶前差
分为
∆n y = ∆(∆n−1 y) = ∆[∆(∆n−2 y)]
⋯⋯ = ∆{∆⋯[∆(∆y)]} = ∆{∆⋯[∆( f (x + ∆x) − f (x)]}
n 阶的向后差分、中心差分的型式类似。
（5-6）
函数的差分与自变量的差分之比，即为函数对自变量的差商。如一阶向前差
二阶差商多取中心式，即
∆2 y ∆x 2
=
f (x + ∆x) − 2 f (x) + (∆x) 2
f (x − ∆x) 。
（5-9） （5-10） （后的二阶差商。 以上是一元函数的差分与差商。多元函数 f(x,y,…)的差分与差商也可以类推。
如一阶向前差商为
应地，上式中的 ∆y 、 ∆x 分别称为函数及自变量的差分， dy //#######为函数对 dx
自变量的差商。 在导数的定义中 ∆x 是以任意方式趋近于零的，因而 ∆x 是可正可负的。在差
分方法中， ∆x 总是取某一小的正数。这样一来，与微分对应的差分可以有 3 种
形式： 向前差分 向后差分 中心差分
高精度傅里叶有限差分法模型正演
高精度傅里叶有限差分法模型正演刘文革 3 ① 贺振华 ①② 黄德济 ① 杜增利 ③( ①成都理工大学信息工程学院 ; ②成都理工大学“油气藏地质与开发工程”国家重点实验室 ; ③西南石油大学)刘文革 ,贺振华 , 黄 德 济 , 杜 增 利 . 高 精 度 傅 里 叶 有 限 差 分 法 模 型 正 演 . 石 油 地 球 物 理 勘 探 , 2007 , 42( 6 ) : 629～633摘要 波动方程数值模拟方法分为有限差分法和频率 —波数域法两类 ,其中有限差分法的计算精度取决于波场外推算子的近似程度 、离散网格间距及差分方程阶数 ,它能适应速度任意横向变化 ,但在大倾角处易出现频 散现象及背景噪声 。
频率 —波数域法算法简单 、精度高 、噪声小 ,能适应任意地层倾角情况 ,但不适于速度场的 任意横向变化 。
文中结合有限差分法和频率 —波数域法的优点 ,应用傅里叶有限差分法 ( F F D) 实现在多域用高 精度延拓算子对模型进行地震记录的数值模拟 ,其波场外推算子由相移项 、折射项 ( 时移项) 和有限差分补偿项 组成 。
对 F F D 法进行了理论与误差分析 ,并用单程声波方程分别进行了层状模型和 SE G / EA GE 盐丘模型的数 值模拟试验 。
数值试验的对比分析表明 , F F D 法适用于速度场横向剧烈变化情形 ,且具有精度高 、无频散 、背景 噪声弱等优点 ,模拟结果反射特征清楚 ,能对复杂地质构造进行准确的地震数值模拟 。
关键词 波动方程 地震数值模拟 有限差分法 频率 —波数域法 傅里叶有限差分法提高微分方程对单平方根算子的逼近程度 ,改善了 成像效果[ 3 ]; C lae r b o u t 在研究时 —空域和频率 —空间域中的有限差分法后 ,认为有限差分法可以适应横向速度的任意变化 ,但是单平方根算子展开引起 了高波数近似 ,使得对大倾角构造的成像误差较大 , 易出现频散现象及背景噪声[ 4 ] 。
时域有限差分法PPT课件
vg
d
dk
c
(1-10)
这种情况下，群速也是与频率无关。
.
8
1.2 数值色散关系(2)
上述过程也可用于一维标量波动方程差分近似的数值色散分析。
设在离散空间点 xi,tn,离散行波解为 u in u x i,tn e j n t k ~ i x ，
式中，k~ 为存在于有限差分网格中的数值正弦波的波数。一般情况 下，不同于连续物理波的波数。正是这种不同导致了数值相速和群 速偏离了精确解。进而导致了数值色散误差。
1.5 数值稳定性（1）
• FDTD计算中每一步都是有误差的，随着时间步进，误 差会不断积累。如果误差的积累不会造成总误差的增 加，就成FDTD法是稳定的，否则成为不稳定的。数值 不稳定性会造成计算结果随时间步进无限增加。
• FDTD法是有条件稳定的，即：时间步必须必须小于一 定值以避免数值不稳定性。
考虑(1.1)的正弦行波解 ux,tejtkx 代入(1-1)得
j2c2jk2 即
k c
上式便是一维标量波动方程的色散关系。
(1-8)
由上式得相速度
vp
k
c
（1-9）
可见，相速与频率无关，称为非色散。非色散意味着对于具有任意
调制的包络或脉冲形状的波传播任意距离后波形保持不变。进一步
由(1-8)可以得到群速关系
正弦函数
ui=sin(nt+)
高斯函数
ui=exp[-(n-n0)2/T2]
阶跃函数
ui= 0
n<n1
= ( n-n1)/(n2-n1) n1<n<n2
=1
n>n2
“硬源”设置简单，但当反射波回到“硬源”位置时， 会引起寄生反射，所以，要在这之前“关”掉源。
2有限差分法及热传导数值计算PPT演示课件
t1
1 a11
(b1
a12t2
a13t3 )
t2
1 a 22
(b2
a 21t1 a 23t3 )
1 t3 a 33 (b3 a 31t1 a 32t2 )
•24
（2）假设一组解（迭代初场），记为： t1(0)、t2(0)并、t代3(0) 入迭代方程求得第一 次解
每次计算t1(1)均、t用2(1)、最t3(1新) 值代入。
(1) 平直边界上的节点
如图所示 边界节点 (m,n) 只能代表半个元体，若边界上有向 该元体传递的热流密度为q ，据能量守恒定律对该元体有：
tm1,n tm,n ytm,n1tm,n x
x
y 2
x
2
tm,n1 tm,n y
Φm,n
2xyyqw
0
Байду номын сангаас
xy tm ,n1 4 2 tm 1 ,n tm ,n 1 tm ,n 1 x2 Φ m ,n 2 x q w
非稳态项 的离散有三种不同的格式。如果将函数在节 点(n，i+1)对点(n，i)作泰勒展开，可有
•30
•31
由式(b)可得在点(n，i)处一阶导数的一种差分表示式 , 的向前差分:
类似地，将t在点(n，i-1)对点(n，i)作泰勒展开，可得 的向后差分的表达式：
如果将t在点(n，i+1)及(n，i-1)处的展开式相加，则可得 一阶导数的中心差分的表达式：
qw
y x
•16
(3) 内部角点
如图所示内部角点代表了 3/4 个元体，在同样的假设条 件下有
tm1,ntm,ny tm,n1tm,nx tm,n1tm,n x
x
有限差分法地震正演模拟程序
有限差分法地震正演模拟程序一．二阶公式推导1．二维的弹性波动方程22222222x x x ZU U U U t x z x z λμμλμρρρ∂∂∂∂++=++∂∂∂∂∂ 22222222xz z z U U U U t z x x zλμμλμρρρ∂∂∂∂++=++∂∂∂∂∂ 2．对方程进行中间差分（1）首先对时间进行二阶差分()()()()()112,,,22222n n n xi k xi k xi kx x x x U U U U t t U t U t t U t t t +--++-+-∂==∂ ()()()()()112,,,22222n n n zi k zi k zi kz z z z U U U U t t U t U t t U t t t +--++-+-∂==∂ （2）对x 方向进行二阶差分()()()()()21,,1,22222n n n xi k xi k xi kx x x x U U U U x x U x U x x U x x x +--++-+-∂==∂ ()()()()()21,,1,22222n n n zi k zi k zi kz z z z U U U U x x U x U x x U x x x +--++-+-∂==∂ （3）对z 方向进行二阶差分()()()()()2,1,,122222n n nxi k xi k xi k x x x x U U U U z z U z U z z U z z z +--++-+-∂==∂ ()()()()()2,1,,122222n n nzi k zi k zi k z z z z U U U U z z U z U z z U z z z +--++-+-∂==∂ （4）对xz 进行差分2,1,11,11,11,11,111,11,11,11,122224n nzi k zi k z zn n n nzi k zi k zi k zi k n n n n zi k zi k zi k zi k U U U U x z x zx z U U U U U U U U x x z z x+-++-++--++-++---⎛⎫-∂∂∂∂⎛⎫== ⎪ ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭-----+==2,1,11,11,11,11,111,11,11,11,122224n nxi k xi k x x n n n n xi k xi k xi k xi k n n n n xi k xi k xi k xi k U U U U x z x z x z U U U U U U U U x x z z x+-++-++--++-++---⎛⎫-∂∂∂∂⎛⎫== ⎪ ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭-----+==（5）把（1）（2）（3）（4）得到的结果带入波动方程3．写出差分方程()()()11,,,1,,1,22,1,,11,11,11,11,12222=2++4n n n n n nxi k xi k xi kxi k xi k xi k nnnnnnnxi k xi k xi k zi k zi k zi k zi k U U U U U U t x U U UU UUUz xz λμρμλμρρ+-+-+-++-++----+-++-+--++()()()11,,,1,,1,22,1,,11,11,11,11,12222=2++4n n n nnnzi k zi k zi kzi k zi k zi knn n nnnnzi k zi k zi k xi k xi k xi k xi k U U U U U U t x U UUU UUUz xz λμρμλμρρ+-+-+-++-++----+-++-+--++即得到()()1,,1,,1,,122112,,,1,11,11,11,1222+=2-++4n n n n n nxi k xi k xi kxi k xi k xi k n n n xi k xi k xi k n n n nzi k zi k zi k zi k U U U U U U x z U U U t U U U U z x λμμρρλμρ+-+-+-++-++---⎛⎫-+-++ ⎪ ⎪ ⎪--++ ⎪ ⎪⎝⎭()()1,,1,,1,,122112,,,1,11,11,11,1222+=2-++4n n n n n nzi k zi k zi k zi k zi k zi k n n n zi k zi k zi k n n n nxi k xi k xi k xi k U U U U U U x z U U U t U U U U z x λμμρρλμρ+-+-+-++-++---⎛⎫-+-++ ⎪ ⎪ ⎪--++ ⎪ ⎪⎝⎭二．MATLAB 程序 clear; clc;Nx=200; Nz=300;fi=30;%%%主频t_step=300;%%%%时间采样点dx=10.0;%空间采样间隔——x 方向 dz=10.0;%空间采样间隔——z 方向 dt=0.001;%时间采样间隔——1mslambda=66*1e9;%砂岩拉梅常数lamdamu=44*1e9;%砂岩拉梅常数murho=2650;%砂岩密度%%%%%%模型扩展%%%%%vp=zeros(Nz+2,Nx+2);%纵波速度vs=zeros(Nz+2,Nx+2);%横波速度c=zeros(Nz+2,Nx+2);%交叉项系数for i=1:Nz+2for j=1:Nx+2vp(i,j)=sqrt((lambda+2*mu)/rho);vs(i,j)=sqrt((mu)/rho);c(i,j)=sqrt((lambda+mu)/rho);endend%%%%%% 震源%%%%%%%t_wavelet=[1:t_step]*dt-1.0/fi;source=((1-2*(pi*fi*t_wavelet).^2).*exp(-(pi*fi*t_wavelet).^2));% 雷克子波amp=sqrt(2.0);% 振幅source_x=floor(Nx/2)+1;% 震源位置——x坐标source_z=2;% 震源位置——z坐标source_amp=zeros(Nz+2,Nx+2);% 震源振幅初始化（所有点处均为0）source_amp(source_z,source_x)=amp;% 震源振幅，在位置source_z,source_x处振幅为amp，其它位置为0 %%%%%%%%%%%%%%%%%%%%%%%%%%%U=zeros(Nz,Nx);% 弹性波x分量V=zeros(Nz,Nx);% 弹性波z分量U0=zeros(Nz+2,Nx+2);% 前一时刻的UU1=zeros(Nz+2,Nx+2);% 当前时刻的UU2=zeros(Nz+2,Nx+2);% 下一时刻的UV0=zeros(Nz+2,Nx+2);% 前一时刻的VV1=zeros(Nz+2,Nx+2);% 当前时刻的VV2=zeros(Nz+2,Nx+2);% 下一时刻的Vrecord_u=zeros(t_step,Nx);% x方向接收到的地震记录——Urecord_v=zeros(t_step,Nx);% x方向接收到的地震记录——V%%%%%% 有限差分计算U V %%%%%%for it=1:t_stepfor x=2:Nx+1for z=2:Nz+1U2(z,x)=2*U1(z,x)-U0(z,x)+(dt*dt)*(vp(z,x)^2*(1.0/(dx*dx))*(U1(z,x+1)-2*U1(z,x) +U1(z,x-1))+vs(z,x)^2*(1.0/(dz*dz))*(U1(z+1,x)-2*U1(z,x)+U1(z-1,x))+c(z,x)^2*(1. 0/(4*dx*dz))*(V1(z+1,x+1)-V1(z+1,x-1)-V1(z-1,x+1)+V1(z-1,x-1)))+source(it)*sour ce_amp(z,x)*dt*dt;V2(z,x)=2*V1(z,x)-V0(z,x)+(dt*dt)*(vs(z,x)^2*(1.0/(dx*dx))*(V1(z,x+1)-2*V1(z,x) +V1(z,x-1))+vp(z,x)^2*(1.0/(dz*dz))*(V1(z+1,x)-2*V1(z,x)+V1(z-1,x))+c(z,x)^2*(1. 0/(4*dx*dz))*(U1(z+1,x+1)-U1(z+1,x-1)-U1(z-1,x+1)+U1(z-1,x-1)));endendU0=U1;U1=U2;V0=V1;V1=V2;record_u(it,:)=U2(2,[2:201]);record_v(it,:)=V2(2,[2:201]);U=U2(2:301,2:201);V=V2(2:301,2:201);end%%%%% 绘制U 和V 的波场图%%%%%%figureimagesc(U(1:300,1:200));title('U');xlabel('x');ylabel('z');figureimagesc(V(1:300,1:200));title('V');xlabel('x');ylabel('z');%%%%% 绘制U 和V 的地震记录%%%%%%figureimagesc(record_u(1:300,1:200));title('U');xlabel('x');ylabel('t_step');figureimagesc(record_v(1:300,1:200));title('V');xlabel('x');ylabel('t_step');U分量波场图V分量波场图U分量地震记录V分量地震记录。