C++实现合成地震记录

我们知道计算合成地震记录的基本原理是,合成地震记录=子波与反射系数的褶积所以需要子波和反射系数.但是用于计算的数据一般是深度域的,要转换到时间域来必须有时深关系.所以.需要的数据:时间/深度关系数据:checkshot或者DT,用于计算反射系数的数据,一般是DT和密度(RHO B).基本步骤:1, 加载数据:如果是斜井的话,加载井斜,计算出SSTVD,设置成Prefered DS(deviation survey);如果有来自VSP或者其他可信渠道的时深关系的话加载进来,叫checkshot,就是时间,深度关系对,用于提供时深关系;加载DT,RHOB曲线;2,数据质量检查:查看checkshot数据覆盖范围,和品质;查看DT,RHOB曲线的品质,如果不好需要用well-edit或者synthetics里带的一些功能进行编辑.DT,RHOB曲线应该是做过Depth match,需要拼接的话是splice好的.3,制作合成地震记录:点击Post,依次选择时深关系,声波曲线,密度曲线(如果没有密度曲线或者品质不好也可以使用经验公式来代替),声波阻抗,反射系数,子波,合成地震记录,地震数据.软件完全是根据原理走的,如果时深关系没选,后续工作无法开展,如果没有DT,密度,就无法生成声波阻抗和反射系数...软件自带有Ricker30经验子波.如果效果不好可以自己提取子波,也可以使用时变子波.4,对比合成地震记录和井旁道实际地震记录,通过bulkshift或者拉伸压缩来调整时深关系.有时需要用c heckshot来校正DT.一般可能先使用Ricker30子波试一下,看看大致情况,如果效果不好,再尝试提取子波.这是一个反复实验的过程.合成地震记录的品质和制作的数据来源的品质有关,对比的好坏和实际地震数据的品质也有关系.总是实际情况总是复杂的.。

利用有限差分方法合成地震记录program seismicimplicit nonereal::lamda1,mu1,r1,dt,dx,dz,vp1,vs1,freq,lamda2,mu2,r2,vp2,vs2real::c11,c12,c21,c22,c31,c32,d1,d2,d3,pireal,dimension(101)::x,zreal,dimension(101,101)::uz1,uz2,uz3,u3,ux1,ux2,ux3integer::i,j,k,n,ccharacter*10::text,nfile1,nfile2,tpwrite(*,*)'input ytpe(mid-reflect,len-direct)'read(*,'(a)') tpif(tp.eq.'mid') thenprint *, 'dx,dz'read(*,*) dx,dzelse if(tp.eq.'len') thenprint *, 'dx,dz'read(*,*) dx,dzend ifvp1=1500.0;vs1=1100.0dt=0.00025freq=10.0r1=300.0lamda1=(vp1**2-2*vs1**2)*r1mu1=vs1**2*r1pi=atan(1.0)*4.00c11=0.5*(lamda1+2.0*mu1)/r1c21=0.5*(lamda1+mu1)/r1c31=mu1/(2.0*r1)d1=(dt/dx)**2d2=dt**2/(4.0*dx*dz)d3=(dt/dz)**2n=500ux1=0.0;uz1=0.0ux2=0.0;uz2=0.0ux3=0.0;uz3=0.0do i=2,100do j=2,100ux2(i,j)=ux1(i,j)+c11*d1*(ux1(i+1,j)+ux1(i-1,j)-2*ux1(i,j)) & +c21*d2*((uz1(i+1,j+1)-uz1(i-1,j+1))-(uz1(i+1,j-1)-uz1(i-1,j-1))) &+c31*d3*(ux1(i,j+1)+ux1(i,j-1)-2*ux1(i,j))uz2(i,j)=uz1(i,j)+c11*d1*(uz1(i,j+1)+uz1(i,j-1)-2*uz1(i,j)) & +c21*d2*((ux1(i+1,j+1)-ux1(i-1,j+1))-(ux1(i+1,j-1)-ux1(i-1,j-1))) &+c31*d3*(uz1(i+1,j)+uz1(i-1,j)-2*uz1(i,j))end doend doux2(51,51)=0.0uz2(51,51)=0.0ux2(50,50)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)uz2(50,50)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)ux2(52,50)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)uz2(52,50)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)ux2(50,52)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)uz2(50,52)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)ux2(52,52)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)uz2(52,52)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)ux2(50,51)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)ux2(52,51)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)uz2(51,50)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)uz2(51,52)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)do k=2,ndo i=2,100do j=2,100ux3(i,j)=-ux1(i,j)+2*ux2(i,j)+2*c11*d1*(ux2(i+1,j)+ux2(i-1,j)-2*ux2(i,j)) & +2*c21*d2*((uz2(i+1,j+1)-uz2(i-1,j+1))-(uz2(i+1,j-1)-uz2(i-1,j-1)))&+2*c31*d3*(ux2(i,j+1)+ux2(i,j-1)-2*ux2(i,j))uz3(i,j)=-uz1(i,j)+2*uz2(i,j)+2*c11*d1*(uz2(i,j+1)+uz2(i,j-1)-2*uz2(i,j)) & +2*c21*d2*((ux2(i+1,j+1)-ux2(i-1,j+1))-(ux2(i+1,j-1)-ux2(i-1,j-1)))&+2*c31*d3*(uz2(i+1,j)+uz2(i-1,j)-2*uz2(i,j))u3(i,j)=sqrt(ux3(i,j)**2+uz3(i,j)**2)end doend doux3(51,51)=0.0uz3(51,51)=0.0ux3(50,50)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)uz3(50,50)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)ux3(52,50)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)uz3(52,50)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)ux3(50,52)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)uz3(50,52)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)ux3(52,52)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)uz3(52,52)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)ux3(50,51)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)ux3(52,51)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)uz3(51,50)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)uz3(51,52)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)write(text,'(i3)') kc=mod(k,10)if(c .eq. 0) thennfile1=text(1:3)//'ux.dat'open(88,file=nfile1)do i=1,101do j=1,101x(i)=(i-1)*dxz(j)=(j-1)*dzwrite(88,'(1x,i3,f10.4,3x,f10.4,3x,f16.4,3x,f16.4,3x,f16.4)')k,x(i),z(j),ux3(i,j),& uz3(i,j),u3(i,j)end doend doclose(88)end ifdo i=2,100do j=2,100ux1(i,j)=ux2(i,j)ux2(i,j)=ux3(i,j)uz1(i,j)=uz2(i,j)uz2(i,j)=uz3(i,j)end doend doend dostopend program seismic实例采用网格数为100*100，x方向和z方向的间隔dx=dz=1，时间间隔dt=0.00025s，模拟地质体的纵波速度vp=1500.0m/s，横波速度为vs=1100.0m/s，震源的频率为freq=10.0 Hz，地质体的密度为r=300.0.震源放在模型区域的中心位置，其振动为爆炸震源。

应用合成地震记录来标定地震层位是地震资料解释中非常重要的手段，也是将地震资料与测井资料相结合的一条纽带。

它最终使抽象的地震数据与实际的地质模型连接起来，为地震资料解释的可靠性提供了依据。

合成记录的精度将直接影响到地震地质层位标定的准确性，因此，提高合成记录的精度就成了地震层位标定的首要问题。

１合成记录的方法原理１．１合成地震记录制作的一般方法一般而言，人工合成地震记录，是利用声波和密度测井资料求取一反射系数序列，再将这一反射系数序列与某一子波反褶积得到结果。

Ｓ（ｔ）　＝　Ｒ（ｔ）　＊　Ｗ（ｔ） （１）式中 Ｓ（ｔ）　——　合成地震记录； Ｒ（ｔ）　——　反射系数序列； Ｗ（ｔ）　——　地震子波。

上式表明，合成记录的好坏与反射系数序列的求取和子波的选择有着密切的关系。

反射系数序列的准确性和精确程度又与测井资料（声波、密度）的采集、处理等过程密切相关；子波的选择，则要考虑子波的长度、相位、频率等诸多因素。

在实际工作中，所得到的结果往往不尽人意［１］，主要表现在：（１）　合成地震记录与井旁地震道附近的地震剖面层位不吻合现象较多，或者说同相轴吻合的时窗长度有限；（２）　合成地震记录与井旁地震道附近的地震剖面能量不吻合现象较多，或者说同相轴“胖瘦”程度吻合有限；（３）　合成地震记录与井旁地震道附近的地震剖面存在一定的时移。

其原因主要在于：①子波受地质条件变化的影响，难以给得恰到好处；②深—时转换存在误差；③褶积模型并不能完全准确地反应地震记录；④实际地震记录存在噪声。

１．２实用优化方法１．２．１校正测井数据首先对测井数据进行校正，对反射系数序列进行非均匀采样［２，３］。

１．２．２选择合适的子波（１）子波的类型。

常用的子波有两类，一是典型子波，如Ｒｉｃｈｅｒ、Ｔｒａｐｅｒｉｏｄ子波等；二是提取子波，提取子波一般有维纳—莱文森混相位子波提取法和自相关子波提取法两种［４，５］。

从剖面提取的实际子波制作的合成记录，虽然其合成地震记录层位精细标定应用研究＊洪余刚 陈景山 代宗仰 李凌峰（西南石油学院资源与环境学院，四川省成都市６１０５００）摘 要：通过对合成记录制作的一般方法进行分析，结合研究区实际地质、地震资料，提出合成记录的制作在层位标定中的实用优化方法，强调了子波的提取方法和子波相位引起的偏差。

一个地震记录道形成的程序/* bed thickness: H=200msampling interval: t=1ms - t=2msvelocity: v=2000m/s - v=5000m/smain frequency: f0<=250Hz*/#include<stdio.h>#include<math.h>#define N 5#define M 40#define pi 3.1415926void main(){FILE *fp;//定义数组V[N]、P[N]、S[N]，分别存储每层P波速度、每层介质密度、每层S波速度floatV[N]={2000,4100,3500,4900,3300,4000},P[N]={2000,2800,2300,2950, 2250,2400},S[N]={1750,4050,2900,4800,2950,3820};//数组Arp[N-1]、R[M+200*(N-1)]，分别表示单独每一层界面处反射P波的反射系数A// R(R[200]、R[400]、R[600]、R[800]、R[1000])为模型中某一层的放射系数//deltaP[N-1]、deltaV[N-1]、deltaS[N-1]，表示相邻两层P波速度、密度、S波速度的差float Arp[N-1],R[M+200*(N-1)],deltaP[N-1],deltaV[N-1],deltaS[N-1];//RW雷克子波float RW[M+200*(N-1)],G[M+200*(N-1)]={0};//入射角angle，主频f0，取样间隔tfloat angle,sum,f0=25,t=0.001;int i,j;//提示输入入射角（0-90）printf("输入入射角(0-90)：");scanf("%f",&angle);//计算层间P波速度、密度、S波速度的差for(i=1;i<N;i++){deltaP[i-1]=P[i]-P[i-1];deltaV[i-1]=V[i]-V[i-1];deltaS[i-1]=S[i]-S[i-1];}//计算反射系数for(i=0;i<M+200*(N-1);i++)//先初始化RR[i]=0;//特殊情形angle=90if(angle==90){for(i=0;i<N-1;i++)Arp[i]=(P[i+1]*V[i+1]-P[i]*V[i])/(P[i+1]*V[i+1]+P[i]*V[i]);for(i=0;i<N-1;i++){if(i==0)sum=1;else{sum=1;for(j=0;j<i;j++)sum*=(1-pow((Arp[j]),2));}R[200*(i+1)-1]=Arp[i]*sum;}}//利用Zoeppritz方程的解，求解反射系数else{angle=angle*pi/180;//switch to radian systemfor(i=0;i<N-1;i++)R[200*(i+1)-1]=(deltaP[i]/P[i]+deltaV[i]/V[i])*1/2+(-deltaP[i]/P[i]+deltaV[i]/V[i]-2*deltaS[i]/S[i])*pow(sin(angle),2)*1/2; }//离散化雷克子波for(i=0;i<M;i++)RW[i]=(1-pow((pi*f0*(i*t)),2))*exp(-pow((pi*f0*(i*t)),2));for(i=M;i<M+200*(N-1);i++)RW[i]=0;//褶积for(i=0;i<(M+200*(N-1));i++){for(j=0;j<(M+200*(N-1));j++)if((i-j)>0&&(i-j)<(M+200*(N-1)))G[i]+=R[j]*RW[i-j];}//将结果输出到"result.txt"if((fp=fopen("result.txt","w+"))==NULL) {printf("open error!\n");}for(i=0;i<M+200*(N-1);i++)fprintf(fp,"%f ",G[i]);。

人工合成地震记录程序设计（一）、人工合成地震记录原理：地震记录上看到的反射波波形是地震子波在地下各反射界面上发生反射时形成的。

反射波的振幅有大有小（决定于界面反射系数的绝对值）、极性有正有负（取决于反射系数的正负）、到达时间有先有后（取决于反射界面的深度）的地震反射子波叠加的结果。

如果地震子波的波形用S (t )表示，地震剖面的反射系数为双程垂直反射时间t 的函数，用R (t )表示，那么反射波地震记录形成的物理过程在数学上就可以用S (t )的R (t )的褶积表示，即某一时刻的反射波地震记录f (t )是:)()()(t R t S t f *=其离散形式为：))(()()(1t m n R t m S t n f Mm ∆-⋅∆=∆∑=如果大地为多层介质，在地面记录长度内可接收的反射波地震记录为：))(()()(11t m n R t m S t n f Mm N n ∆-⋅∆=∆∑∑== 式中，n 为合成地震记录的采样序号，n =1,2,3...N ；N 为合成一道地震记录的采样点数；m =1,2,3...M ，为离散子波的采样点数；△t 为采样间隔。

这种褶积模型将地震波的实际传播过程进行了简化：1、在合成地震记录的过程中没有考虑大地的吸收作用，所有薄层的反射波都与地震子波的形式相同，只是振幅和符号不同。

2、假设地震波垂直入射到界面上，并原路径返回。

3、假设地层横向是均匀的，在深度（纵向）方向上假设密度为常数，只是速度发生变化。

4、不考虑地震波在传播过程中的透射损失。

（二）、人工合成地震记录的方法1、 反射系数序列在有速度测井资料的情况下，可以用速度曲线代替波阻抗曲线，计算反射系数序列。

在没有速度资料的情况下，可根据干扰波调查剖面分析的结果设计地质模型。

如设计的地质模型如图a 所示，图中H 为层厚度，V 为层速度，根据下式计算反射系数： 11)(--+-=N N N N N V V V V H R 式中H 为反射界面的深度，N 为反射层序号，随深度变化的反射系数序列如图b 所示。

地震合成记录
地震合成记录是指通过模拟地震波的传播过程，生成一组地震波形数据，以便研究地震活动的特征以及对建筑物和基础设施的影响。

这些合成记录可以用于地震工程、地震学和地质学等领域的研究和应用。

在进行地震合成记录时，首先需要确定地震事件的源参数，包括震级、震源深度、断层类型以及断层面和滑动方向等信息。

接着，利用地震波传播理论和数值模拟方法，计算出地震波在地下的传播路径和幅度衰减规律。

最后，将计算得到的地震波形数据与已知的地震记录进行比较和校正，以确保模拟结果的准确性。

地震合成记录的应用非常广泛。

在地震工程中，合成记录可以用于评估建筑物和结构在地震中的响应，预测地震对建筑物的破坏程度，以及指导抗震设计和加固工作。

在地震学领域，合成记录可以用于研究地震波的传播特性，深入了解地球内部的结构和物理特性，以及预测地震的发生概率和强度。

在地质学领域，合成记录可以用于研究地震活动对地层的影响，揭示地壳运动和地震活动的地质机制。

随着计算机技术的不断进步，地震合成记录的方法和技术也在不断发展。

近年来，基于高性能计算和大规模并行计算的地震波传播模拟已经取得了显著的进展，使得合成记录的精度和时空分辨率得到了提高。

此外，随着地震数据库的不断积累和更新，研究人员可以更准确地选择合适的地震事件和地震记录进行合成，提高合成结果的可靠性和适用性。

地震合成记录在地震研究和工程实践中扮演着重要的角色。

通过合成记录的研究和分析，可以更好地理解地震的力学过程，提高地震预测和防灾减灾能力，为社会的安全和可持续发展提供支持。