《反演与层析成像》实验

合集下载
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

本科生实验报告
实验课程地球物理反演与层析成像
学院名称地球物理学院
专业名称地球物理学
学生姓名
学生学号
指导教师
实验地点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++)
{
ray[12*i+j].begin=1.5+3.0*i;
ray[12*i+j].end=1.5+3.0*j;
ray[12*i+j].slope=(ray[12*i+j].begin-ray[12*i+j].end)/45;
ray[12*i+j].length=sqrt(45*45+pow(fabs(ray[12*i+j].begin-ray[12*i+j ].end),2));
}
}
fp=fopen("射线.txt","w");
for(i=0;i<144;i++)
{
fprintf(fp,"%d %f %f %f %f\n",i+1,ray[i].begin,ray[i].end,ray[i].sl ope,ray[i].length);
}
fp1=fopen("慢度.txt","w");
for(i=0;i<12;i++)
{
for(j=0;j<9;j++)
fprintf(fp1,"%f ",xc[i][j]);
fprintf(fp1,"\n");
}
for(i=0;i<144;i++)
for(j=0;j<12;j++)
for(y=0;y<9;y++)ray[i].transform[j][y]=0;
for(i=0;i<144;i++)
{
for(j=0;j<100;j++)
{
ray[i].point[j].x=0;
ray[i].point[j].y=0;
}
}
for(i=0;i<144;i++)
{
ray[i].time=0;
}
fp4=fopen("y.txt","w");
for(i=0;i<12;i++) {
for(j=0;j<12;j++)
{
for(q=0;q<100;q++)
{
d[q].x=100;
d[q].y=0;
}
for(m=0;m<10;m++)
{
d[m].x=m*5.0;
d[m].y=-d[m].x*ray[i*12+j].slope+ray[i*12+j].begin;
}
n=10;
e=ray[i*12+j].n;
if(e>0)
{
for(y=1;y<=e;y++)
{
d[n].y=ray[i*12+j].begin+3.0*y-1.5;
d[n].x=(3.0*y-1.5)/fabs(ray[12*i+j].slope);
n++;
}}
if(e<0)
{
for(y=-1;y>=e;y--)
{
d[n].y=ray[i*12+j].begin+3.0*y+1.5;
d[n].x=-(3.0*y+1.5)/fabs(ray[12*i+j].slope);
n++;
}}
b=0;
while(d[b].x!=100)
{
b++;
}
for(c=0;c<b;c++)
for(k=0;k<b-1;k++)
if(d[k].x>d[k+1].x)
{
temp1=d[k].x;
d[k].x=d[k+1].x;
d[k+1].x=temp1;
temp2=d[k].y;
d[k].y=d[k+1].y;
d[k+1].y=temp2;
}
b=0;
while(d[b].x!=100)
{
ray[12*i+j].point[b].x=d[b].x;
ray[12*i+j].point[b].y=d[b].y;
fprintf(fp4,"%d\t%f,%f
\n",12*i+j+1,ray[12*i+j].point[b].x,ray[12*i+j].point[b].y);
b++;
}
for(a=1;a<b;a++)
{
mm=sqrt(pow((ray[12*i+j].point[a-1].x-ray[12*i+j].point[a].x),2)+pow((ray[1 2*i+j].point[a-1].y-ray[12*i+j].point[a].y),2)) ;
if(ray[12*i+j].n>=0)
ray[12*i+j].transform[(int)(ray[12*i+j].point[a-1].y/3)][(int)(ray[
12*i+j].point[a-1].x/5)]=mm;
if(ray[12*i+j].n<0)
ray[12*i+j].transform[(int)(ray[12*i+j].point[a].y/3)][(int)(ray[12
*i+j].point[a-1].x/5)]=mm;
}
}
}
for(i=0;i<144;i++)
{
for(j=0;j<12;j++)
{
for(y=0;y<9;y++)
{
ray[i].time+=ray[i].transform[j][y]*xc[j][y];
}
}
}
fp3=fopen("矩阵A.txt","w");
for(i=0;i<144;i++)
{
for(j=0;j<12;j++)
{
for(y=0;y<9;y++)
fprintf(fp3,"%f ",ray[i].transform[j][y]);
fprintf(fp3,"\n");
}
fprintf(fp3,"\n");
}
fp5=fopen("矩阵A2.txt","w");
for(i=0;i<144;i++)
{
for(j=0;j<12;j++)
{
for(y=0;y<9;y++)
fprintf(fp5,"%f ",ray[i].transform[j][y]);
}
fprintf(fp5,"\n");
}
fp6=fopen("走时.txt","w");
for(i=0;i<144;i++)
{
fprintf(fp6,"%f\n",ray[i].time);
}
}
三、结果图
实验六射线追踪-向后处理
一、实验原理
在向后进行反向追踪时,也要考虑来自各个方向的射线。

找出来自把该接收点围起来的各个边界射线中,到达该接收点旅行时最小的射线,确定此射线来自的边界以及边界上的插值点,把该插值点作为下一个接收点,以此类推,至到接收点为炮点,就完成了一个接收点到炮点的射线追踪。

在此过程中,要记录每条射线经过每个格子的长度。

实验七~八反演方法
一、实验原理
计算理论走时
利用每条射线经过每个格子的长度以及经过格子的慢度,便可求得每条线的理论走时。

上述的改进虽然可以解决原LTI算得的节点旅行时不是最小的问题,却降低了计算效率。

在一些复杂的模型里,它的效率会是原LTI的1/5甚至更低。

所以为了兼顾精度和效率的要求,可以设定循环次数。

另外,在比较粗糙的初始模型初始迭代可以用原LTI来进行射线追踪,到后面反演得到的模型精度比较高之后,再用改进LTI。

此外,在实际工程当中,每炮和它的接收点对应的射线只会覆盖整个区域很小的一部分,所以为了提高效率,在对每一炮进行射线追踪时,可以只计算这些射线覆盖的区域,而不必计算整个区域。

这个区域可以由炮点最左边的接收点和最右边的接收点来确定,考虑到有回波的可能,此区域可适当地扩大。

下面是计算循环次数不相同得到的射线追踪图。

原模型是在坐标为(1500m ,300m)到(4000m,300m),以及从(0,700m)到(1500m,700m)有一水平的高速反射界面。

在(1500m,300m)到(1500m,700m)有一坚直的相同速度的高速反射界面,水平反射界面以上是从地表初速度800m/s以每增深10m,速度增大30m/s。

炮点在原点,图上所标的S的位置。

每个接点相隔100m,第一个接收点在炮点右侧200m处。

(a) 循环0次
(b) 循环10次
(c) 循环20次
(d) 循环到没有节点的旅行时减小为止
图4
从上图可以看出,算节点旅行时时循环次数不同,追踪出来的射线也不同。

随着循环次数的增多,出现了明显的回波现象。

表1给出每条射线在循环次数为0,10,20以及循环到没有节点的旅行时减小为止计算出来的走时。

(有“\”的地方表示走时没有变化)
表1
同上表可以看出,Rec13, Rec14,Rec15对应的走时经多次循环计算后有明显的减小,循环到没有节点的旅行时减小为止算得的走时分别比没有循环计算得到的走时小27.16ms,81.46ms,9.04ms 。

而这些接收点接收到的地震波是遇到高速反射界面后反射回来的。

2.4 模型修正
用LTI 进行射线追踪得到每个接收点在初始模型上的理论走时后,就利用这个理论走时与实际拾取走时的差值,对初始模型进行修正。

以下是本文用SIRT 算法修正模型的处理步骤:
1,假设初始模型的第j 个网格中的慢度为j S ,第k 个炮点对应的接收点有
k N 个;
其中,j=1,2,…,M(M 为网格总数); k=1,2,…,NUM (NUM 为炮点总数) 2,利用前面的射线追踪方法,得到该炮点的每个接收点的理论走时n P
1
M
n nj j
j P a S ==∑, n=1,2…,k N (2.1)
式中nj a 由射线追踪得到的第n 条射线在第j 个网格内射线长度
3,求出该炮点的每个接收点实际拾取走时n T 与理论走时
n P 之差
n t ∆
n n n t T P ∆=- , n=1,2…,k N (2.2)
4, 设第n 条射线在第j 个网格内的慢度修正值
nj c ,则
1
M
n j
n j n
j a c t
==∆∑, n=1,2…,k N (2.3
设修正值nj c
正比于第j 个网格内射线通过的路径nj a 与该射线长n R 之比,即:
nj nj n
n
a c R α= ,n=1,2…,k N (2.4)
式中n α是第n 条射线的比例常数,n R 是第n 条射线的全长
1
M
n n j j R a ==∑, n=1,2,…,k
N
(2.5)
将(2.5)(2.4)式代入(2.3)式,并整理简化,可得
2
1
()
nj
nj n
M
nj
j a c t a
==∆∑ , n=1,2,…,k N (2.6)
1
k
N j j nj n c c c ==+∑, (2.7)
其中,j c 为第j 个网格的累记修正量,在算第一个炮点之前初始化为0 5, 重复步骤2,3,4,计算第k+1个炮点,直到计算完NUM 个炮点。

6,求每个网格的平均修正值,设计算完NUM 个炮点后,经过第j 个网格的
射线总条数为j Y ,则每个网格的平均修正值为
j j j
c c Y =
,n=1,2,…,k N (2.8)
7,用平均修正值对第 j 个 网格的慢度值j S 进行修正,即:
j j j S S c =+, (2.9)
在修正后,j S 值需受下列物理条件的约束,即
min max j S S S ≤≤, (2.10)
若j S <min S ,则取j S =min S ; 若j S >max S ,则取j S =max S ;
min S , max S 为介质慢度的范围,不同的介质,慢度值的范围不同。

直到这里,便完成了一次迭代。

从上面的计算步聚可以看出,SIRT 不是用一条射线的修正值nj c 来对慢
度j S 进行修正的,而是将所有射线得到的修正值保存下来,在本轮对射线迭
代结束后求所有射线在网格内的修正值 的平均值来对慢度j S 来修正的。

为了提高解的可靠性,对穿过射线总条数少于3条的网格的慢度不作处理,而是用以下的插值方法得到。

插值处理
因为射线是有限的,不可能所有的网格里都会有射线经过,没有射线经过的网格,则该网格的慢度就得不到修正。

没有修正值的网格则可以通过周围的有修正值的网格慢度插值得到。

本文采用的插值方法是反距离平方法,即把未知点与已知点的距离作为权重因子,未知点与已知点的距离越近,其权重越大,反之越小,权重由距离平方的反比给出,表达式为:
211
[][]m i i i
S S d d ==∑∑m
2i=1i (2.11)
式中,S 为未知点,即需要插值得到的网格慢度;i S 为已知点,即周围已修正过的第i 个网格慢度;i d 是未知点到第i 个已知点的距离,m 为参与插值
的已知点的个数。

参与插值的已知点分布范围不应太大,否则加长搜索时间,降低计算效率,而且离未知点较远的已知点的值对未知点的值影响已经很小了。

可以设定一个半径,与未知点的距离在这个半径范围内的已知点才参与插值。

插值处理不是任意情况都适用。

如果已知点比未知点多,而且未知点是离散地分布在已知点中间,则这些未知点插值得到插值得到的值可靠性就高。

如果未知点比已知点多,或者未知点是连续分布在一大块区域内,则这些未知点插值得到插值得到的值可靠性就低。

对修正后的模型进行平滑处理
为了使更新的速度场变化不过于激烈,更接近于实际当中的速度场分布特点,以及保证反演的稳定性,有必要对更新的速度场作平滑处理。

本文的平滑处理方法如下:每个网格的慢度等于以该网格为中心、大小为一定值的矩形区域内所有慢度的平均值。

矩形区域不宜太大,太大会降低精度,减慢收敛速度。

矩形区域太小,平滑的效果不佳。

本文采用的是3×3个网格的矩形区域。

2.7 假设第q 次迭代得到的j S 用()
q j S 来表示,对求得的j S 值,用下式判断其收敛程度
()(1)
||q q j j
S S ε--< (2.12) 如果条件(2.11)式成立,则认为j S 值达到预定的收敛要求,否则以此时的j S 值作为初始值再做下一轮近迭代,直到满足(2.11)式的条件时停止。

反演近地表速度结构
本文有幸得到中国石油勘探研究院提供的中国西部地区西秋的一条二维测
线的地震资料。

该测线位于山地地区,测线长约为33000m,共计348炮,每炮各有200~500个道。

为了验证算法的有效性,本文以实际资料的地形设置理论模型作为反演结果的参考标准,如图5所示。

炮与道的位置,井深均按照实际资料来设置。

本文的主要任务是在该理论模型做正演,得到理论的拾取走时,通过反演把浅层的速度在2500m/s 以下的低降速带(图中黑线以上部分)反演出来。

本文计算时采用的初始模型如图6所示,该模型从地表开始,以1800m/s 的初速度,每增深10米,速度增大15m/s 的梯度往下递增。

网格大小为32m ×10m 。

图5 反演的理论模型
图6 反演的初始模型
炮检距的影响
本文着重解决的是近地表速度结构的反演问题。

野外观测系统提供了足够大的炮检距,但反演近表速度结构时并不是取的炮检距越在越好。

虽然炮检距越大反演的深度越大,射线数越多,但由于计算过程中的平均效应,深层的误差会带到表层来,从而使反演的表层速度的精度降代。

所以如果只为反演较浅的表层速度结构,则只需小炮检距的初至信息即可。

图7(a)和图7(b)分别为最大炮检距取2000m和3500m反演迭代5次的层析结果图。

图7(b)的反演深度比图7(a)的大,但精度却不如图7(a)。

图8是最大炮检距2000m和3500m层析结果地表以下50m的速度切片。

与理论模型的相对比,迭代相同的次数,最大炮检距2000m的更接近于理论模型。

图9显示的是这两个最大炮检距目标函数值随迭代次娄的变化,可以看出炮检距小的收敛得更快一下。

这里得出的结论跟刘玉柱等人的是一样的。

但并不是炮检距越小越好。

如果覆盖近表地层的射线严重不足,也会影响反演精度。

最大炮检距可以根据炮间距,在初始模型上正演时射线走的深度以及想要的反演深度来确定。

图7(a) 最大炮检距2000m层析结果
图7 (b) 最大炮检距3500m层析结果
图8 最大炮检距2000m(蓝线)、最大炮检距3500m(蓝线)层析结果与理论模型对比
速度切片(地表以下50米)
图9最大炮检距2000m (蓝线)和3500m (绿线)的层析收敛曲线
改进SIRT 算法
受到炮检距问题的启发,本文还对SIRT 算法进行了改进和发展。

SIRT 算法求网格的修正值时,是对通过该网格的所有射线对应的修正值取平均来得到的,也就认为的所有射线对应的修正值的地位是一样的。

从SIRT 算法求修正值来看,对于同一条射线穿过的网格,走时误差n t ∆按射线穿过每个网格内长度来分配,从而得每个网格的修正值。

可走时误差有可能只由射线穿过的部分网格的慢度误差造成的,如果按SIRT 那样算,没有慢度误差的网格也给修正了。

虽然这样能减小n t ∆,但修正得到的速度场却是与实际偏离了的。

所以穿得越深或走的路径越长的射线对其经过的网格的修正值就越不可靠。

基于这点考虑,本文认为每条射线对网格的修正值的可靠性是不一样的:射线越长,经过的格子越多,可靠性就越低;射线越短,经过的格子数越少,可靠性就越高。

所以本文在求出所有射线对网格的修正值后,乘上一个权重系数,射线经过的总格子数越少,权重越大。

设第n 条射线经过的总格子数为n W ,权重由射线经过的总格子数的反比给出时,,则(2.7)式改写为:
1
j j nj n
c c c W =+⋅ ,n=1,2,…, k
N (3.1)
若权重由射线经过的总格子数的平方的倒数给出时,则(2.7)式改写为:
2
1
j j nj n
c c c W =+⋅
上式的n 也可以只取nj c >0对应的n ;设nj c >0对应的n W ,i=1,2,…, j Y ,
权重由射线经过的总格子数的反比给出时,(2.8)式改写为: 11j j
j Y i n c c W ==∑ (3.2)
若权重由射线经过的总格子数的平方的倒数给出时,则(2.8)式改写为: 211j j
j Y i n c c W ==∑
对慢度进行修正时,仍然按照(2.9)式计算。

图10(a)和图10(b)为最大炮检距均为2000m ,迭代次数均为5次,加权系数分别由射线经过的总格子数的倒数给出(加权1)和由总格子数平方的倒数给出(加权2)时的层析结果。

与图7(a)的相比,加权处理后的反演精度更高,而加权2反演出来的近地表速度结构更接近于理论模型。

图11显示的是层析结果地表以下50m 的速度切片,从该图也可以看出,加权处理比没作加权处理的精度高,而加权2比加权1的精度高。

图12显示的是没有加权,加权1,加权2的层析收敛曲线,由图可以看出加权处理的比没加权处理的收敛得更快一些,而加权2又比加权1的要快。

图10(a )加权系数由射线经过的总格子数的倒数给出时层析结果(加权1)
图10 (b) 加权系数由射线经过的总格子数平方的倒数给出时层析结果(加权2)
图11 没有加权(紫红线)、加权1(绿线)、加权2(蓝线)层析结果与理论模型(红线)
对比
地表以下50m速度切片
图12 没有加权(蓝线),加权1(绿线),加权2(红线)的层析收敛曲线
实际地震资料处理
为了验证以上的分析结果,本文还对上述二维测线实际的地震走时数据进行了反演计算。

初始模型如图6所示,网格大小仍为32m×10m,最大炮检距取2000m,采用加权2处理。

图13(a)是该测线的近地表速度结构的层析结果。

图13(b)显示了该测线部分炮点的拾取初至(红线)与在层析结果上计算的初至(蓝线)的对比。

从该图可以看出理论初至与拾取初至总体上吻合得比较好。

图14是应用层析结果进行静校正的效果图。

图14(a)为未应用层析静校正剖面图,图14(b)为应用层析静校正剖面图,两图对比来看,应用了层析静校正后,剖面图中有改进的地方,但也有变差的地方。

说明层析结果的精确度仍有所欠缺,需要进一步提高。

分析其原因可能有以下几个:一是地震资料中的初至拾取是有误差的,而在反演计算时却当作是无误差的;二是在二维测线中所有的炮和道的排列并不是严格地在一条直线上的,可在计算时却看作所有的炮和道都在一条直线上;三是用LTI模拟计算地震波在地下从炮点传到接收点的初至时间与实际的存在偏差,数学模拟是把实际情况理想化了,而实际情况往往会更复杂;四是用SIRT 进行反演求解时,虽然能满足使走时残差减到很小,但并不一定满足解的唯一性。

五是求静校正量的方法有待进一步探讨。

图13 二维实际测线最终层析结果
图13(b )二维测线部分炮拾取初至(红线)与理论计算初至(蓝线)对比
(a )未应用层析静校正剖面 (b )应用层析静校正剖面 图14 层析静校正效果图
本次实习体会
本文讨论了用LTI-SIRT 反演近地表速度结构的方法和步骤。

LTI 是一种精度高和速度快的射线追踪方法,但存在不能追踪逆向传播的射线路径。

经本文的改进后,克服了这一缺点,但同时也使计算效率降低。

但随着计算机的速度越来越快,这也将不成为问题。

此外,讨论了炮检距对反演近地表速度结构的影响,通过理论模型的验算,发现炮检距太大时,虽然反演深度增大,但把深层的误差带到浅层来,从而影响反演浅层速度的精度。

与此同时,本文还对SIRT 进行了改进,认为射线越短对网格速度的修正可靠性高,越长则可靠性越低,从而在对速度进行修正时加上权系数,而不是简单地取平均。

通过理论模型的验算,改进后的方法可以提高层析收敛速度和反演精度。

从理论模型的计算中可以看出,用LTI-SIRT
反演近地表速度结构得到的结果趋势是正确的,但在细节上还存在偏改进
变差
差。

本文还用实际地震初至拾取数据进行层析反演,并利用层析结果作静校正来验证层析结果的好坏。

从应用层析静校正剖面图来看,有改进的地方,也有变差的地方。

造成此结果的可能原因是多方面的。

山地的近地表速度构造往往是很复杂的,要很好地反演出此构造得依靠野外观测系统误差的减小,先验信息的利用,算法适应性高,参数的适当选择等。

相关文档
最新文档