《反演与层析成像》实验

本科生实验报告

实验课程地球物理反演与层析成像

学院名称地球物理学院

专业名称地球物理学

学生姓名

学生学号

指导教师

实验地点5417

实验成绩

二〇-四年十二月二〇一五年一月

实验一射线追踪-建立观测系统

一、实验原理

1.1 初至拾取

初至拾取也就是对初至波的到达时间的记录，这是进行层析反演计算的基础数据。数据的好坏直接关系到层析成像的效果好坏。如果是进行实际资料计算时，初至拾取则为野外采集所得数据；若为理论模型的计算，则是在理论模型上进行正演计算所得的初至波旅行时。

1.2 建立初始模型

初始模型的建立，就是按照实际模型的地表起伏给定一假设的速度场。为方便起见，这个速度场可以是从地表开始以某一初始速度以相同的梯度往下递增。这一初始速度可以同先验知识给出，又或者从初至拾取中获得。炮－道比较接近时，可以把地震波看作直线传播，已知炮检距及初至时间就可以算出地表速度。初始模型的深度一定要够大，能够让射线自由地传播，不会出现遇到模型底部强迫反射的现象。

1.2 射线追踪

射线追踪的准确与否是影响层析成像的关键。鉴于日本科学家Aszkawar提出的LTI算法的高效率和高精度，本文采用该算法进行旅行时计算及射线追踪。LTI算法基于Fermat原理,即地震波沿着一条传播时间最短的路径进行传播。

该方法把模型离散成均匀的正方形单元，旅行时和射线路径的确定只与单元边界上的点有关。假设单元边界上任一点的旅行时可由该边界上相邻两个离散点的旅行时线性插值得到。如图1所示为一匀速正方形单元，单元边界平行于坐标轴，A,B是二个旅行时己知的点，要求射线穿过A,B边界到达D点的最小旅行时及射线路径。设射线从C点通过，C点旅行时可用线性内差公式由A,B二点的值表示。D点的旅行时为C点旅行时与波在C, D间直线传播时间之和。然后根据Fermat原理，就可求出C点的位置（即射线路径）和D点的最小旅行时。该方法也分成向前和向后处理，向前处理只计算各单元边界上的节点的旅行时；向后处理根据Fermat原理追踪射线路径．确定的射线路径不是单元边界上离散节点的连线，而是穿过单元边界上正好满足最小旅行时条件的那一点的线。这样，在一个单元内射线路径为直线，同时边界上的折射角随入射角连续变化。所以，

该方法具有较高的精度和计算速度。

图1匀速正方形单元旅行时线性插值几何关系

A、B 为已知点，C是内插点，D是待求点

二、程序源代码

#include"stdio.h"

#include"math.h"

struct data

{ double begin;

double end;

double slope;

double length;

double transform[12][9];

double transform2[108];

int n;

double time;

struct point

{

double x;

double y;

}point[100];

}ray[144];

struct ww

{

double x;

double y;

}d[100];

实验二~五射线追踪-向前处理

一、实验原理

考虑到要反演地下速度结构，对纵向的分辨率要求要高，所以网格纵向长度要比横向长度小。网格取得越小，分辨率越高，但计算量和数据存储量同时也会加大，所以要综合考虑计算效率以及内存因素。本文在对模型离散化时，采用矩形网格，计算节点只取网格线的交点。本文采用的网格模型如图2所示：

图2

LTI存在的问题是，不能追踪到逆向传播的射线路径。这是因为，原LTI的计算方向只考虑了初至波正向传播的情况，而没有考虑逆向传播的情况，当模型速度变化复杂时，就不能追踪逆向传播的初至射线。此外，原LTI算得的节点旅行时并不一定都是最小的。原LTI算每列未知的节点旅行时时，算的节点先后顺序不一样，就可能算出不一样的值，而且无法确定先算哪个节点后算哪个节点。因为根据最小走时原理，应该是先算地震波先到达的节点也即旅行时小的节点，后算地震波晚到达的节点。可原LTI算法一开始无法确知地震波先到达哪个节点后到达哪个节点。为了解决这些问题，本文对原LTI作出了改进，具体步骤为：1，向前处理：

（1）按原LTI得到各节点的旅行时；

（2）设置一个记录节点旅行时是否有变化的标志数组，数组元素对应每个节点，初始值都为1；

（3）除了炮点所在网格边上的节点，其余的节点再算一次。此时插值的边界包括把该节点围起来的各个边界，如图2－3所示，如果边界两端点

的标志都为0，或者两端点的旅行时都比该节点的大，则不再计算此边

界。如果算到比原来小的，就用新算出来的值替代原来的。

图3

（4）比较此次计算得到的每个节点的旅行时和上一次计算得到的节点旅行时，如果不相等，则置与该节点的标志为1，如果相等,则置为0；

（5）重复步骤（3），（4），至到没有节点的旅行时不再变小为止。此时便能得到节点的最小旅行时。

二、程序源代码

main()

{ FILE *fp,*fp0,*fp1,*fp2,*fp3,*fp4,*fp5,*fp6;

int i,j,y,m,b,e,n,q,c,k,a;

double mm,temp1,temp2;

double x[108],xc[12][9];

for(i=0;i<108;i++)x[i]=1.0/3000.0;

x[20]=x[29]=1.0/5000.0;

x[77]=x[78]=1.0/2000.0;

for(i=0;i<12;i++)

for(j=0;j<9;j++)xc[i][j]=x[9*i+j];

fp0=fopen("shijimoxing.txt","w");

for(i=0;i<12;i++)

{

for(j=0;j<9;j++)

fprintf(fp0,"%f",xc[i][j]);

fprintf(fp0,"\n");

}

for(i=0;i<12;i++)

{

for(j=0;j<12;j++)

{