《反演与层析成像》实验 - 图文

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

一、 实验原理

考虑到要反演地下速度结构，对纵向的分辨率要求要高，所以网格纵向长度要比横向长度小。网格取得越小，分辨率越高，但计算量和数据存储量同时也会加大，所以要综合考虑计算效率以及内存因素。本文在对模型离散化时，采用矩形网格，计算节点只取网格线的交点。本文采用的网格模型如图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(\ for(i=0;i<12;i++) { }

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

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

fprintf(fp0,\

fprintf(fp0,\

{

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\ for(i=0;i<144;i++) {

fprintf(fp,\ope,ray[i].length); }

fp1=fopen(\慢度.txt\ for(i=0;i<12;i++) {

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

fprintf(fp1,\}

fprintf(fp1,\ }

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(i=0;i<144;i++) { }

ray[i].time=0; for(j=0;j<100;j++) {

ray[i].point[j].x=0; ray[i].point[j].y=0; }

fp4=fopen(\

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

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

{

for(q=0;q<100;q++) { }

for(m=0;m<10;m++) { } n=10;

d[m].x=m*5.0;

d[m].y=-d[m].x*ray[i*12+j].slope+ray[i*12+j].begin; d[q].x=100; d[q].y=0;

e=ray[i*12+j].n;

if(e>0) { {

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++;

for(y=1;y<=e;y++)

}} if(e<0) { {

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); for(y=-1;y>=e;y--)

n++; }} b=0;

while(d[b].x!=100)