地震数值模拟实验报告

本科生实验报告

实验课程数值模型模拟

学院名称地球物理学院

专业名称勘查技术与工程

学生姓名ZRY

学生学号

指导教师

实验地点624

实验成绩

二〇一五年4月二〇一五年5月

成都理工大学

《地震数值模拟》实验报告

实验二叠加地震记录的相移波动模拟实方程正演验

摘要

利用C语言编制地质模型的相移波动方程正演模拟，改变绕射点位置、速度，再做正演模拟。

关键字：地震模型；正演记录

1.1实验目的

掌握各向同性介质任意构造、水平层状速度结构地质模型的相移波动方程正演模拟基本理论、实现方法与程序编制，由正演记录初步分析地震信号的分辨率。

1.2实验内容

1、基本要求：

（1）点绕射构造和水平层状速度模型（参数如图1 所示）的正演数值模拟；

1）削波的正演；

2）无削波的震正演；

（2）计算中点和两个边界的信号位置，分析实验结果的正确性；

（3）做同样模型的褶积模型数值模拟，对比分析分析两者的异同。

（4）改变绕射点位置、速度，再做正演模拟。

2、较高要求：

（1）使用雷克子波做爆炸源，对三个不同的主频：25hz、50hz 和75hz 分别做点绕射

模型的正演模拟；

（2）设计复杂反射构造模型，再做正演模拟。

1.3实验原理

1、地震波传播的波动方程

设（x,z）为空间坐标，t 为时间，地震波传播速度为v(x,z)，则二位介质中任意位置、任意时刻的地震波场为p(z,x,t)：压缩波——纵波。则二维各向同性均匀介质中地震波传播

的遵循声波方程为

（）

2、傅里叶变换的微分性质

p(t)与其傅里叶变换的P(ω)的关系：

则有时间微分性质

ω 为频率，ω=2π/T，T 为周期。同理有空间微分性质：

k 为频率，k=2π/λ，λ为波长。

3、地震波传播的相移外推公式

令速度v 不随x 变化，只随z 变化，则利用傅里叶变换微分性质（3）和（4）式，把波动方程(1)式变换到频率-波数域，得：

或：

令：

则（5）式的解为：

包括上行波和下行波两项。正演模拟取上行波：

若和间隔为△Z ，速度v(z) 在此间隔内不随Z 变的常数,(7)式实现波场从到的延拓，即：

在深度Zj+1 开始向上延拓到Zj，若延拓深度为零，即：∆Z= Z j+1-Z j=0，则

对于任意深度Z j+1 到Z j 的延拓，可得正演模拟中地震波的传播方程（延拓公式）

4、初始条件和边界条件

按照爆炸界面理论，反射界面震源在t=0 时刻同时起爆，此时刻的波场就是震源。根据不同情况，可直接使用反射系数脉冲或子波作震源。如果直接使用反射系数作震源脉冲，则初始条件可表示为：

其他

对时间t和空间x做二维傅立叶变换，则得频率-波数域的初始波场

边界条件：

其他，其他，其他

其他参数都是在范围内定义的。

5、边界处理

（1）边界反射问题

把实际无穷空间区域中求解波场的问题化为有穷区域求解时，左右两边使用零边界条件。物理上假设探区距和两个端点很远，在两个端点上收到的反射波很弱。但是，上述条件在实际中不能成立，造成零边界条件反而成为绝对阻止波通过的强反射面。在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。

（2）边界强反射的处理

镶边法、削波法、吸收边界都能有效消除边界强反射。

削波法就是在波场延拓过程中，没延拓一次，在其两侧均匀衰减到零，从而消除边界强反射的影响。假设横向总长度为NX，以两边Lx 道吸波为例，有以下吸波公式：

其他

6、数字化

根据数字信号处理的采样定理，把连续的信号变为计算机能处理的数字信号，使相移法正演模拟得以实现。

频域抽样定理：一个频谱受限信号f(t)，如果时间只占据-t m—+t m 的范围，若在频域以不大于1/2t m 频率间隔∆f≤1/2t m 对信号f(t)的频谱F(ω)采样，则抽样到的离散信号F1(ω)可以唯一表示原信号。

时域抽样定理：一个时间受限信号f(t)，如果频谱只占据-ωm—+ωm 的范围，则信号f(t)可以用等间隔的抽样值唯一表示出来，而时间∆t 抽样间隔必须不大于1/2f m，ωm =2π f m，∆t≤1/2f m。

1.4实验过程

//#include "PsFrwrdMdlParameter.h"

#include

#include

#include

#include

#include

int kkfft(float pr[], float pi[], int n, int l);

int Absorb(int Labs);

int Po2Judge(int);

int Rflct();

int Vlcty();

int WvFld0();

int PsFrwd();

int exp_ikzDz(float eikzdz[],int Ix, float Vc,int Iw,float Dw,float Dkx);

int PhaseShift();//相移

int Frqcy2Time();

#define Nx 128 //Trace Number

#define Nt 256 //Record Number