声波方程有限差分正演

声波方程有限差分正演
声波方程有限差分正演

题目:使用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)

对分解后的声波方程进行离散,可得到:

1

12211,-1,,,122[]N n n n n m i m j i m j xi j xi j m t V V c P P h +

-+---=?=+-∑ 1

1

221

1,1,,,122

[]N n n n n m i j m i j m zi j zi j m t V V c P P h +-++---=?=+-∑ 111121

2222,,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),使用二阶中心差商公式:

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 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 变换:

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 z w t t c h

??+?=?) (5) 公式(5)右端必须满足下列条件:

2

2222sin sin 220(x z

k x k z t c h

??+≤?≤)1 取x k 和z k 最大值,即=x x k π?,z =k z π?,则有:22

220t 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 ,一个波长内的空间采样点数 100m=1010h λ=

=,基本满足网格色散条件。 程序运行输出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 = newfloat * [n];

for (int i = 0; i

{

A[i] = newfloat [n];

}

float *x = newfloat [n];

float *b = newfloat [n];

//

for (int i = 0; i

{

for (int j = 0; j

{

A[i][j] = pow( 2.0*(j+1)-1, 2.0*(i+1)-1 );

}

}

for (int i = 0; i

{

if (i< 1)

b[i] = 1.0;

else

b[i]=0;

}

//

Gauss(A,n,b,x);

for (int i = 1; i<= n; i++)

{

c1[i] = x[i-1];

}

for (int i = 0; i

{

delete[] A[i];

}

delete A;

delete x;

delete b;

}

主函数:

void main()

{

FILE *fp;

char name[1000];

float dt, h;

float Ts, f0;

int NX, NZ, NT, s_x, s_z, nPoints;

fp=fopen("2D_Parameters.txt", "r");

fscanf(fp, "%[^\n]\n", name);

fscanf(fp, "%f\n", &dt); //Time Interval

fscanf(fp, "%[^\n]\n", name);

fscanf(fp, "%d\n", &NX); //X Grid Dimension

fscanf(fp, "%[^\n]\n", name);

fscanf(fp, "%d\n", &NZ); //Z Grid Dimension

fscanf(fp, "%[^\n]\n", name);

fscanf(fp,"%d\n",&s_x); //Source X

fscanf(fp, "%[^\n]\n", name);

fscanf(fp, "%d\n", &s_z); //Source Z

fscanf(fp, "%[^\n]\n", name);

fscanf(fp,"%f\n",&h); //Space Interval

fscanf(fp, "%[^\n]\n", name);

fscanf(fp, "%f\n", &Ts); //Total Time

fscanf(fp, "%[^\n]\n", name);

fscanf(fp, "%f\n", &f0); //Dominant Frequency

fscanf(fp, "%[^\n]\n", name);

fscanf(fp, "%d\n", &nPoints); //Accuracy of FD method

int model;

fscanf(fp, "%[^\n]\n", name);

fscanf(fp, "%d\n", &model); //model

fclose(fp);

NT = (int)((Ts + 1e-6)/dt);

///////////////////////////////////////////////////////////////// //The process of opening the array

float **U1, **U2;

float **Txx1, **Txx2;

float **Tzz1, **Tzz2;

float **Vp;

U1 = Create2DActivityArray(NZ,NX);

U2 = Create2DActivityArray(NZ,NX);

Txx1 = Create2DActivityArray(NZ,NX);

Txx2 = Create2DActivityArray(NZ,NX);

Tzz1 = Create2DActivityArray(NZ,NX);

Tzz2 = Create2DActivityArray(NZ,NX);

Vp = Create2DActivityArray(NZ,NX);

// Set Model

Velocity_Model(model, Vp, NX, NZ);

float *c1;

c1 = (float*)malloc( sizeof(float) * (nPoints+1) );

Cal_1D_FdCoff(c1, nPoints);// Compute FD Coefficient

// Print Coefficient

for (int i = 1; i<= nPoints; i++)

{

printf("%e\n", c1[i]);

}

///////////////////////////////////////////////////////////////// printf("***********Forward Calculation starts!************\n");

float sum1, sum3;

for (int k = 0; k

{

for (int j = 0; j

{

for (int i = 0; i

{

sum1 = 0.0;

sum3 = 0.0;

for (int m = 0; m

{

// U/X

if ( i-m-1 < 0 )

{

sum1 += (float)c1[m+1] * Txx1[j][i+m];

}

elseif ( i+m>NX-1 )

{

sum1 -= (float)c1[m+1] * Txx1[j][i-m-1];

}

else

{

sum1 += (float)c1[m+1] * (Txx1[j][i+m] -

Txx1[j][i-m-1]);

}

// U/Z

if ( j-m< 0 )

{

sum3 += (float)c1[m+1] * Tzz1[j+m+1][i];

}

elseif ( j+m+1 >NZ-1 )

{

sum3 -= (float)c1[m+1] * Tzz1[j-m][i];

}

else

{

sum3 += (float)c1[m+1] * (Tzz1[j+m+1][i] -

Tzz1[j-m][i]);

}

}

//

U2[j][i] = -Vp[j][i]*Vp[j][i] * (dt/h) * (sum1 + sum3) + U1[j][i];

}//for i

}//for j

///////////////////////////////////////////////////////////////// for (int j = 0; j

{

for (int i = 0; i

{

sum1 = 0.0;

sum3 = 0.0;

for (int m = 0; m

{

// V/X

if ( i-m< 0 )

{

sum1 += (float)c1[m+1] * U2[j][i+m+1];

}

elseif ( i+m+1 >NX-1 )

{

sum1 -= (float)c1[m+1] * U2[j][i-m];

}

else

{

sum1 += (float)c1[m+1] * (U2[j][i+m+1] -

U2[j][i-m]);

}

// V/Z

if ( j-m-1 < 0 )

{

sum3 += (float)c1[m+1] * U2[j+m][i];

}

elseif ( j+m>NZ-1 )

{

sum3 -= (float)c1[m+1] * U2[j-m-1][i];

}

else

{

sum3 += (float)c1[m+1] * (U2[j+m][i] -

U2[j-m-1][i]);

}

}

//

Txx2[j][i] = -(dt/h) * sum1 + Txx1[j][i];

Tzz2[j][i] = -(dt/h) * sum3 + Tzz1[j][i];

//添加震源项

if ( i==s_x&&j==s_z )

{

Txx2[j][i]+=f(k*dt, f0);

}

}//for i

}//for j

////////////////////////////////////////////////////////////

////Tranfer the value.////

for (int j = 0; j

for (int i = 0; i

{

U1[j][i] = U2[j][i];

Txx1[j][i] = Txx2[j][i];

Tzz1[j][i] = Tzz2[j][i];

}

}

////////////////////////////////////////////////////////////

if (k> 0)

{

Fwrite_Snapshot_U(U2, NX, NZ, k, 50);

}

printf("The %d ms Forward file has finished!\n",k);

//Finish the circulation of K.

}

Release2DActivityArray(U1,NZ);

Release2DActivityArray(U2,NZ);

Release2DActivityArray(Txx1,NZ);

Release2DActivityArray(Txx2,NZ);

Release2DActivityArray(Tzz1,NZ);

Release2DActivityArray(Tzz2,NZ);

Release2DActivityArray(Vp,NZ);

}

(完整版)大连理工大学高等数值分析抛物型方程有限差分法

抛物型方程有限差分法 1. 简单差分法 考虑一维模型热传导方程 (1.1) )(22x f x u a t u +??=??,T t ≤<0 其中a 为常数。)(x f 是给定的连续函数。(1.1)的定解问题分两类: 第一,初值问题(Cauchy 问题):求足够光滑的函数()t x u ,,满足方程(1.1)和初始条件: (1.2) ()()x x u ?=0,, ∞<<∞-x 第二,初边值问题(也称混合问题):求足够光滑的函数()t x u ,,满足方程(1.1)和初始条件: ()13.1 ()()x x u ?=0,, l x l <<- 及边值条件 ()23.1 ()()0,,0==t l u t u , T t ≤≤0 假定()x f 和()x ?在相应的区域光滑,并且于()0,0,()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。

现在考虑边值问题(1.1),(1.3)的差分逼近 取 N l h = 为空间步长,M T = τ为时间步长,其中N ,M 是 自然数, jh x x j ==, ()N j ,,1,0Λ=; τ k y y k ==, ()M k ,,1,0Λ= 将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。其中 ()j i y x ,表 示网格节点; h G 表示网格内点(位于开矩形G 中的网格节点)的集合; h G 表示位于闭矩形G 中的网格节点的集合; h Γ表示h G -h G 网格边界点的集合。 k j u 表示定义在网点()k i t x ,处的待求近似解,N j ≤≤0,M k ≤≤0。 注意到在节点()k i t x ,处的微商和差商之间的下列关系 ((,)k j k j u u x t t t ????≡ ? ????): ()() ()ττ O t u t x u t x u k j k j k j +??? ????=-+,,1 ()() ()2112,,ττ O t u t x u t x u k j k j k j +??? ????=--+ ()()()h O x u h t x u t x u k j k j k j +??? ????=-+,,1 ()() ()h O x u h t x u t x u k j k j k j +??? ????=--,,1 ()() ()2112,,h O x u h t x u t x u k j k j k j +??? ????=--+ ()()() ()2 222 11,,2,h O x u h t x u t x u t x u k j k j k j k j +???? ????=+--+ 可得到以下几种最简差分格式

叠加地震记录的相移波动方程正演模拟数值模拟实验共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道吸波为例,有以下吸波公式: 四、实验内容

【毕业设计(论文)】二维热传导方程有限差分法的MATLAB实现

第1章前言 1.1问题背景 在史策教授的《一维热传导方程有限差分法的MATLAB实现》和曹刚教授的《一维偏微分方程的基本解》中,对偏微分方程的解得MATLAB实现问题进行过研究,但只停留在一维中,而实际中二维和三维的应用更加广泛。诸如粒子扩散或神经细胞的动作电位。也可以作为某些金融现象的模型,诸如布莱克-斯科尔斯模型与Ornstein-uhlenbeck过程。热方程及其非线性的推广形式也被应用与影响分析。 在科学和技术发展过程中,科学的理论和科学的实验一直是两种重要的科学方法和手段。虽然这两种科学方法都有十分重要的作用,但是一些研究对象往往由于他们的特性(例如太大或太小,太快或太慢)不能精确的用理论描述或用实验手段来实现。自从计算机出现和发展以来,模拟那些不容易观察到的现象,得到实际应用所需要的数值结果,解释各种现象的规律和基本性质。 科学计算在各门自然科学和技术科学与工程科学中其越来越大的作用,在很多重要领域中成为不可缺少的重要工具。而科学与工程计算中最重要的内容就是求解科学研究和工程技术中出现的各种各样的偏微分方程或方程组。 解偏微分方程已经成为科学与工程计算的核心内容,包括一些大型的计算和很多已经成为常规的计算。为什么它在当代能发挥这样大的作用呢?第一是计算机本身有了很大的发展;第二是数值求解方程的计算法有了很大的发展,这两者对人们计算能力的发展都是十分重要的。 1.2问题现状 近三十年来,解偏微分方程的理论和方法有了很大的发展,而且在各个学科技术的领域中应用也愈来愈广泛,在我国,偏微分方程数值解法作为一门课程,不但在计算数学专业,而且也在其他理工科专业的研究生的大学生中开设。同时,求解热传导方程的数值算法也取得巨大进展,特别是有限差分法方面,此算法的特点是在内边界处设计不同于整体的格式,将全局的隐式计算化为局部的分段隐式计算。而且精度上更好。 目前,在欧美各国MATLAB的使用十分普及。在大学的数学、工程和科学系科,MATLAB

抛物形扩散方程的有限差分法及数值实例

偏微分方程数值解 所在学院:数学与统计学院 课题名称:抛物形扩散方程的有限差分法及数值实例学生姓名:向聘

抛物形扩散方程的有限差分法及数值实例 1.1抛物型扩散方程 抛物型偏微分方程是一类重要的偏微分方程。考虑一维热传导方程: 22(),0u u a f x t T t x ??=+<≤?? (1.1.1) 其中a 是常数,()f x 是给定的连续函数。按照初边值条件的不同给法,可将(1.1.1)的定解分为两类: 第一,初值问题(Cauchy 问题):求足够光滑的函数()t x u ,,满足方程(1.1.1)和初始条件: ()()x x u ?=0,, ∞<<∞-x (1.1.2) 第二,初边值问题(也称混合问题):求足够光滑的函数()t x u ,,满足方程(1.1.1)和初始条件: ()()x x u ?=0,, 0x l << (1.1.3) 及边值条件 ()()0,,0==t l u t u , T t ≤≤0 (1.1.4) 假定()x f 和()x ?在相应的区域光滑,并且于()0,0,()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。 1.2抛物线扩散方程的求解 下面考虑如下热传导方程 22()(0.)(,)0(,0)()u u a f x t x u t u L t u x x ????=+????? ==??=??? (1.2.1) 其中,0x l <<,T t ≤≤0,a (常数)是扩散系数。 取N l h = 为空间步长,M T =τ为时间步长,其中N ,M 是自然数,用两族

平行直线jh x x j ==, ()N j ,,1,0 =和k t t k τ ==, ()M k ,,1,0 =将矩形域 G {}T t l x ≤≤≤≤=0;0分割成矩形网格。其中 (),j k x t 表示网格节点;h G 表示 网格内点(位于开矩形G 中的网格节点)的集合;h G 表示位于闭矩形G 中的网格节点的集合;h Γ表示h G -h G 网格边界点的集合。 k j u 表示定义在网点(),j k x t 处的待求近似解,N j ≤≤0,M k ≤≤0。 现在对方程进行差分近似: (一) 向前差分格式 =-+τ k j k j u u 111 2 2(())k k k j j j j j j u u u a f f f x h +--++= (1.2.2) ()j j j x u ??==0, k u 0=k N u =0 (1.2.3) 计算后得: 111(12)k k k k j j j j j u ru r u ru f τ++-=+-++ (1.2.4) 其中,2 a r h τ = ,1,,1,0-=N j ,1,,1,0-=M k 。 显然,这是一个四点显示格式,每一层各个节点上的值是通过一个方程组求解到的。方程组如下: 1000 121011000 232121000 3432310001121(12)(12)(12)(12)N N N N N u ru r u ru f u ru r u ru f u ru r u ru f u ru r u ru f ττττ----?=+-++?=+-++??=+-++? ???=+-++? (1.2.5) 若记 () T k N k k k u u u 1 21,,,-= u ,()()()()T N x x x 121,,,-=???? ,()()()()T N x f x f x f 121,,,-=τττ f 则显格式(1.2.4)可写成向量形式 10 ,0,1,,1 k k k M φ +?=+=-?=? u Au f u (1.2.6) 其中

【文献综述】热传导方程差分格式的收敛性和稳定性

文献综述 信息与计算科学 热传导方程差分格式的收敛性和稳定性在实际研究物理问题过程中, 往往能给出问题相应的数学表达式, 但是由于实际物理问题的复杂性, 它的解却一般不容易求出. 由此计算物理应运而生, 计算物理是以计算机为工具, 应用数学的方法解决物理问题的一门应用性学科, 是物理、数学和计算机三者结合的交叉性学科. 它产生于二战期间美国对核武器的研究, 伴随着计算机的发展而发展. 计算物理的目的不仅仅是计算, 而是要通过计算来解释和发现新的物理规律. 这一点它与传统的实验物理和理论物理并无差别, 所不同的只是使用的工具和方法. 计算物理早已与实验物理和理论物理形成三足鼎立之势, 甚至有人提出它将成为现代物理大厦的“栋梁”. 在一个物理问题中一个数值解往往比一个式子更直观, 更有价值. 在实际求解方程时, 除了一些特殊的情况下可以方便地求得其精确解外, 在一般情况下, 当方程或定解条件具有比较复杂的形式, 或求解区域具有比较复杂的形状时, 往往求不到, 或不易求到其精确解. 这就需要我们去寻找方程的近似解, 特别是数值近似解, 简称数值解. 这里主要研究的是热传导方程. 有限差分法是微分方程和积分微分方程数值解的方法. 其基本思想是把连续的定解区域用有限个离散点构成的网格来代替, 这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似, 积分用积分和来近似, 于是原微分方程和定解条件就近似地代之以代数方程组, 即有限差分方程组, 解此方程组就可以得到原问题在离散点上的近似解. 然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解. 热传导的差分法是求解热传导方程的重要方法之一. 对于差分格式的的求解, 我们首先要关注差分格式的收敛性和稳定性. 对于一个微分方程建立的各种差分格式, 为了有实用意义, 一个基本要求是它们能够任意逼近微分方程, 即相容性要求. 一个差分格式是否有用, 就要看差分方程的精确解能否任意逼近微分方程的解, 即收敛性的概念. 此外, 还有一个重要的概念必须考虑, 即差分格式的稳定性. 因为差分格式的计

声波方程数值模拟实验报告

声波方程数值模拟实验报告 一.基础理论知识 需要的已知条件包括: 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年代,美国哥伦比亚大学在郭宗汾

基于三阶Adams格式的求解声波方程的多步算法

创新项目论文 一种基于三阶Adams 格式的求解声波方程的多步算法 China University of Mining & Technology-Beijing

摘要 一个准确、高效、低数值频散的正演算法能够提高反演精度、加快反演收敛速度,因此研究地震波场正演模拟技术具有重要意义。区别于传统的空间离散方法,利用空间插值, 用网格点处的函数值及其梯度共同逼近空间高阶偏导数的方法称为近似解析离散化方法。声波方程通过变换,并采用近似解析离散化方法进行空间离散,从而转变成为一个半离散化的常微分方程组,再利用三阶显式Adams格式进行时间推进,求解半离散化的常微分方程组,从而得到了一个新的求解声波方程的有限差分方法(AD-STEM)。对AD-STEM进行了理论误差和数值误差分析、计算效率比较和数值波场模拟。研究表明,与传统方法AD-LWC比较,AD-STEM方法数值精度更高,数值频散更低,更高效,且与解析解匹配更好。AD-STEM方法能够通过压制数值频散而提高计算效率。在无可见数值频散的条件下,AD-STEM的计算速度是AD-LWC的1.88倍,而存储量只有其72%,更适合在粗网格下进行大规模地震波场数值模拟。 关键词:近似解析离散化方法;三阶Adams格式;数值频散;有限差分

目录 1 绪论 (1) 1.1选题背景和研究意义 1.2粘弹性介质国内外研究现状 1.3有限差分国内外研究现状 1.4本文主要研究内容 2 粘弹性介质的基本模型 (6) 3方法介绍....................................................................................................................... 错误!未定义书签。 3.1 Stereo-modeling方法简介 (10) 3.2 Lax-Wendroff correction方法简介 ...................................................................... 错误!未定义书签。 4 粘弹性介质中的波场数值模拟..................................................................................... 错误!未定义书签。 4.1 波场快照 (11) 4.2 波形图.................................................................................................................. 错误!未定义书签。 4.3 SEG模型的地表地震记录 (14) 5 结论 (18) 6 参考文献 (20)

2007射线追踪与波动方程正演模拟方法对比研究

47 科技资讯  科技资讯 SCIENCE & TECHNOLOGY INFORMATION2007 NO.12 SCIENCE & TECHNOLOGY INFORMATION 工 业 技 术 地震正演模拟作为反演解释的反过程,是验证解释成果的有效手段,进行必要可靠的正演模拟可以有效的监控反演解释。地震学一般可以分为几何地震学和物理地震学,在几何地震学中进行的正演模拟方法就是我们通常所说的射线追踪法,射线追踪法是在合成记录时用地震子波和界面或地质体的反射系数进行反褶积运算,即。运算的最大特点是说明了地震波传播的运动学特征。而在物理地震学中应用波动方程法合成的地震记录是通过求解波动方程的数值解来模拟地震波场的。在波动方程合成的地震记录中不单保持了地震波传播运动学特征,还说明了地震波传播的动力学特征。本文将分别用射线追踪和波动方程的方法合成地震记录。 1 基于射线追踪的合成地震响应 射线追踪法的主要理论基础是,在高频近 射线追踪与波动方程正演模拟方法对比研究 王志美 畅永刚 (长江大学油气资源与勘探技术教育部重点实验室 湖北荆州 434023) 摘 要:地震学一般可以分为几何地震学和物理地震学,几何地震学中进行正演模拟方法就是射线追踪法,射追踪法是在合成记录时用地震子波和界面或地质体的反射系数进行反褶积运算,即。运算的最大特点就是说明了地震波传播的运动学特征。而在物理地震学中的波动方程法合成的地震记录是通过求解波动方程的数值解来模拟地震波场。在波动方程合成的地震记录中不单保持了地震波传播 运动学特征外,还说明了地震波传播的动力学特征。本文将分别用射线追踪和波动方程的方法合成地震记录。关键词:射线追踪 波动方程 正演模拟 中图分类号:P315文献标识码: A 文章编号:1672-3791(2007)04(c)-0047-00 图1 射线追踪正演模拟(1) 图 2 逐段迭长示意图 图 3 射线追踪正演模拟(2) 图4 波动方程正演模拟结果 似条件下,地震波的主能量沿射线轨迹传播。基于这种认识,运用惠更斯原理和费马原理来重建射线路径,并利用程函方程来计算射线的旅行时。在旅行时计算中应用有限差分等方法,以获得快速的解。射线法的主要优点是概念明确,显示直观,运算方便,适应性强;其缺陷是应用有一定限制条件,计算结果在一定程度上是近似的,对于复杂构造进行两点三维射线追踪往往比较麻烦。为了计算波沿射线的旅行时和波的传播路径,叙述如下。 如图1所示,首先给出连接S(激发点)和R(接收点)之间的初始射线路径射线的振幅变化,首先必须知道地震波在实际地层中传播的射线路径。 由于地震波在整条路径上满足同一个射线参数,因此射线路径上任意连续三点也将满足同一个参数,而三点间的射线表现形式为Snell定律。按照Snell 定律,可导出一个求 取中间点的一阶近似公式。当前后两点位于界面两边时,中间点为透射点,所求路径为透射路径;当前后两点位于界面的同一边时,中间点为反射点,所求路径为反射路径。为此,可以从任一端点出发,连续地选取三点,通过一阶近似公式进行逐段迭代取中间点,利用新求出的点代替原来的点,然后以一点的跨跃作为步长,顺序地逐段迭代下去,直到另一端点。这样,新计算出的中间点和两个端点就构成了一次迭代射线路径,如图2中所示。如果整条射线路径上校正量的范数之和满足一定的精度要求,则认为射线追踪过程结束,否则从追踪出的射线路径开始,继续重复上述过程,直到满足精度要求为止。最后一次追踪到的中间点和两个端点,构成整条射线路径。图3基于多层倾斜界面模型通过射线追踪正演模拟地震响应。从模拟结果可以直观的看出基于几何地震学的原理正演模拟结果只能反映地震波的几何传播路径。在实际的工程设计中通过正演模拟可以在地表确定地下观测范围,节约设备提高工程效率,但不能反映 物理地震学中的地震属性,例如振幅,频率和相位等。更不能反映地震波的动力学特征。 2 波动方程的合成地震响应 2.1 波动方程的建立 非均匀介质的声波方程:  (1) (2) 可由对连续介质方程(1)式的两端对时间求导,并利用欧拉方程推得:  (3) 其中:P是波数,V是质点振动的速度向量,ρ是密度,c是波速,ρ和c是随着空间参数χ和z变化的,这里ρ给定为常数,只有c 是地质模型的控制参数。χ和Z分别是在地面水平距离和深度。这样(3)式就可以变为:  (4) 其中:c=ν (χ,z);(4)式即是所求的弹性波动方程。 2.2 数值计算及稳定性 求解弹性波动方程的方法有多种,付立叶变换法是对弹性波动方程的波场进行付立叶变换,优点是运算速度快。克希霍夫积分法是基于均匀模型,利用格林函数公式计算曲面积分,求出空间波场值,但这种方法不能适应

一维导热方程 有限差分法 matlab实现

第五次作业(前三题写在作业纸上) 一、用有限差分方法求解一维非定常热传导方程,初始条件和边界条件见说明.pdf 文件,热扩散系数α=const , 22T T t x α??=?? 1. 用Tylaor 展开法推导出FTCS 格式的差分方程 2. 讨论该方程的相容性和稳定性,并说明稳定性要求对求解差分方程的影响。 3. 说明该方程的类型和定解条件,如何在程序中实现这些定解条件。 4. 编写M 文件求解上述方程,并用适当的文字对程序做出说明。(部分由网络搜索得到,添加,修改后得到。) function rechuandaopde %以下所用数据,除了t 的范围我根据题目要求取到了20000,其余均从pdf 中得来 a=0.00001;%a 的取值 xspan=[0 1];%x 的取值范围 tspan=[0 20000];%t 的取值范围 ngrid=[100 10];%分割的份数,前面的是t 轴的,后面的是x 轴的 f=@(x)0;%初值 g1=@(t)100;%边界条件一 g2=@(t)100;%边界条件二 [T,x,t]=pdesolution(a,f,g1,g2,xspan,tspan,ngrid);%计算所调用的函数 [x,t]=meshgrid(x,t); mesh(x,t,T);%画图,并且把坐标轴名称改为x ,t ,T xlabel('x') ylabel('t') zlabel('T') T%输出温度矩阵 dt=tspan(2)/ngrid(1);%t 步长 h3000=3000/dt;

h9000=9000/dt; h15000=15000/dt;%3000,9000,15000下,温度分别在T矩阵的哪些行T3000=T(h3000,:) T9000=T(h9000,:) T15000=T(h15000,:)%输出三个时间下的温度分布 %不再对三个时间下的温度-长度曲线画图,其图像就是三维图的截面 %稳定性讨论,傅里叶级数法 dx=xspan(2)/ngrid(2);%x步长 sta=4*a*dt/(dx^2)*(sin(pi/2))^2; if sta>0,sta<2 fprintf('\n%s\n','有稳定性') else fprintf('\n%s\n','没有稳定性') error end %真实值计算 [xe,te,Te]=truesolution(a,f,g1,g2,xspan,tspan,ngrid); [xe,te]=meshgrid(xe,te); mesh(xe,te,Te);%画图,并且把坐标轴名称改为xe,te,Te xlabel('xe') ylabel('te') zlabel('Te') Te%输出温度矩阵 %误差计算 jmax=1/dx+1;%网格点数 [rms]=wuchajisuan(T,Te,jmax) rms%输出误差

声波方程有限差分正演

题目:使用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) 对分解后的声波方程进行离散,可得到: 1 12211,-1,,,122[]N n n n n m i m j i m j xi j xi j m t V V c P P h + -+---=?=+-∑ 1 1 221 1,1,,,122 []N n n n n m i j m i j m zi j zi j m t V V c P P h +-++---=?=+-∑ 111121 2222,,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),使用二阶中心差商公式: 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 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 +--

有限差分法求解抛物型方程说明

有限差分法求解抛物型方程 偏微分方程只是在一些特殊情况下,才能求得定解问题解的解析式,对比较复杂的问题要找到解的解析表达式是困难的,因此需采用数值方法来求解.有限差分法是一种发展较早且比较成熟的数值求解方法,只适用于几何形状规则的结构化网格.它在微分方程中用差商代替偏导数,得到相应的差分方程,通过解差分方程得到微分方程解的近似值.本章主要介绍有限差分法的基本思想,并给出一些具体的数值实例. §1 差分方法的基本思想 有限差分法把偏微分方程的求解区域划分为有限个网格节点组成的网格,主要采用Taylor 级数展开等方法,在每个网格节点上用有限差分近似公式代替方程中的导数,从而建立以网格节点上的函数值为未知数的代数方程组. 有限差分格式,从格式的精度来划分,有一阶格式、二阶格式和高阶格式.从差分的空间形式来考虑,可分为中心格式和逆风格式.考虑时间因子的影响,差分格式还可以分为显格式、隐格式和显隐交替格式等.目前常见的差分格式,主要是上述几种格式的组合,不同的组合构成不同的差分格式. 泰勒级数展开法对有限差分格式的分类和公式的建立起着十分重要的作用.下面采用泰勒展开式导出一个自变量系统的若干有限差分表达式. 首先考虑单变量函数()u x ,如图1把区域x 离散为一批结点,记 0()(), =0,1,2,i i u x u x ih u i =+= 图1 单变量函数离散化 函数()u x 在点i x 处的泰勒展开式为 23 ()()()()()2!3! i i i i i u x u x u x h u x u x h h h ''''''+=++ ++ (1) 或 23 ()()()()()2!3! i i i i i u x u x u x h u x u x h h h ''''''-=-+ -+ (2) 式(1)和(2)重新整理可得 2()()()()()2!3! i i i i i u x h u x u x u x u x h h h '''''+-'= --- (3)

二维频率域声波方程正演模拟

Open Journal of Natural Science 自然科学, 2020, 8(4), 258-263 Published Online July 2020 in Hans. https://www.360docs.net/doc/3e12393984.html,/journal/ojns https://https://www.360docs.net/doc/3e12393984.html,/10.12677/ojns.2020.84034 2D Acoustic Wave Equation Forward Modeling in the Frequency Domain Kun Han, Xiangchun Wang* School of Geophysics and Information Technology, China University of Geosciences (Beijing), Beijing Received: Jun. 23rd, 2020; accepted: Jul. 6th, 2020; published: Jul. 13th, 2020 Abstract Forward modeling in frequency domain plays an important role in the numerical simulation of seismic waves. Compared with time domain forward modeling, frequency domain forward mod-eling has many advantages, such as suitable multi shot parallel operation, no time dispersion, flexible frequency band selection and small error. The coefficient matrix of different frequencies is relatively independent in the frequency domain forward modeling, which is suitable for the acce-leration of parallel computing and greatly improves the computing efficiency. In this paper, for the optimal 9-point difference scheme of frequency domain acoustic equation, the implicit expression and sparse matrix solution are studied, and the seismic wave field is simulated forward. The ac-curacy and validity of the method are verified by model calculation. Keywords Frequency Domain, Forward Modeling, Acoustic Equation, Parallel Computing 二维频率域声波方程正演模拟 韩坤,王祥春* 中国地质大学(北京),地球物理与信息技术学院,北京 收稿日期:2020年6月23日;录用日期:2020年7月6日;发布日期:2020年7月13日 摘要 频率域正演在地震波数值模拟中占有十分重要的地位。相比于时间域正演,频率域正演具有适合多炮并*通讯作者。

基于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)

FDMOD–声波方程有限差分正演模拟二维

FDMOD –声波方程有限差分正演模拟(二维) 格式: fdmod wfile nx= nz= tmax= xs= zs= [optional parameters] 必需的参数: wfile 波场输出文件(包含每个时间步的波场值wave[nx][nz])nx= x采样点个数(第二维) nz= z采样点个数(第一维) xs= 炮点x坐标 dxs= 炮点x坐标间隔 zs= 炮点z坐标 dzs= 炮点z坐标间隔 ns= 炮点个数 tmax= 最大记录时间 可选参数: nt=1+tmax/dt 时间采样点数(dt决定结果的稳定度) mt=1 波场输出时间切片的时间步长间隔 dx=1.0 x采样间隔 fx=0.0 x起始值 dz=1.0 z采样间隔 fz=0.0 z起始值 fmax = vmin/(10.0*h) 震源子波的最高频率 fpeak=0.5*fmax 雷克子波的峰值频率 dfile= 密度输入文件(包含密度值d[nx][nz]) vsx= 垂直测线的x坐标 hsz= 水平测线的z坐标 rsx= 水平测线的起始检波器x坐标 rlen= 水平测线长度 rivl= 水平测线检波器采样间隔 vsfile= 垂直测线的输出文件data[nz][nt] hsfile= 水平测线的输出文件data[nx][nt] ssfile= 震源点检波器的输出文件data[nt] verbose=0 =1 显示输出信息=2 更多输出信息 abs=1,1,1,1 模型的顶,底,左,右使用吸收边界条件 =0,1,1,1 顶部使用自由边界条件

毕设论文--粘声波正演模拟研究

本科毕业设计(论文)题目:粘声波正演模拟方法研究 学生姓名:xxx 学号:xxx 专业班级:xxx 指导教师:xxx 2015年 6月20日

粘声波正演模拟方法研究 摘要 地球上介质的黏滞性会引起大地的吸收效应,它会影响波场所有的频率成分,尤其对于高频的影响最大,导致地震分辨率降低。黏滞吸收作用会影响地震波波形、频带、振幅等因素。一个高效的粘声波正演模拟方法,可以考虑到由于实际介质造成的地震波的吸收衰减作用。可以更加准确模拟地震波在非完全弹性实际地层中的传播,在这里,本文通过编程建立不同的粘声波方程数值模拟模型跟正常的声波方程数值模拟模型进行对比分析,从而了解粘声波正演模拟方法的优越性。 关键词:粘声波;正演模拟;有限差分;

Study on the forward modeling of viscoelastic acoustic waves Abstract The absorption effect is mainly caused by the viscosity of the earth media itself.The viscous stagnation can affect all the frequency components of the wave field.And the effect of the high frequency components is bigger,which leads to the decrease of seismic resolution.The absorption of the absorption has a great influence on the wave, frequency and amplitude of the seismic wave.. A highly effective viscoelastic forward modeling method can take into account the absorption and attenuation of seismic waves by real media.. Accurate simulation of the propagation of seismic waves in the actual strata of the imperfect elasticity. Here. In this paper, the program, establish different visco acoustic wave equation numerical simulation model with normal acoustic wave equation numerical simulation model for comparative analysis, to understand the visco acoustic forward modeling method of superiority. Keywords:Viscoelastic acoustic wave;Viscoelastic acoustic wave;Finite difference;

声波波动方程正演模拟程序总结

声波波动方程正演模拟程序 程序介绍: 第一部分:加载震源,此处选用雷克子波当作震源。 编写震源程序后,我将输出的数据复制,然后我用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阶精度

相关文档
最新文档