地震子波波形显示及一维地震合成记录
- 格式:docx
- 大小:527.92 KB
- 文档页数:22
人工合成地震记录程序设计(一)、人工合成地震记录原理:地震记录上看到的反射波波形是地震子波在地下各反射界面上发生反射时形成的。
反射波的振幅有大有小(决定于界面反射系数的绝对值)、极性有正有负(取决于反射系数的正负)、到达时间有先有后(取决于反射界面的深度)的地震反射子波叠加的结果。
如果地震子波的波形用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 M m ∆-⋅∆=∆∑=如果大地为多层介质,在地面记录长度内可接收的反射波地震记录为:))(()()(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 所示。
地震正演模拟
1. 原理
在地震记录上看到的波形是地震子波叠加的结果,从地下许多反射界面发生反射时形成的地震子波,振幅大小决定于反射界面反射系数的绝对值,极性的正负决定于反射系数的正负,到达时间的先后取决于界面深度和覆盖层的波速。
若地震子波波形用S(t)表示,反射系数是双程垂直反射旅行时t 的函数,用R(t)表示,地震记录f(t)形成的物理过程在数学上就可表示为:
f(t)=S(t)* R(t)=∫S (t )R (t −τ)dτT
其中反射系数R 取R=ρ2v2−ρ1v1ρ2v2+ρ1v1,设地层密度为均匀的,且令ρ2=4ρ1,则反射系数只与地层间速度有关。
地震子波取雷克子波S(t)=[1-2*(pi*fp*t )2]exp[-(pi*fp*t )2],其中主频取fp=20。
2. 模型
本次模拟中采用的是300*80的网格,道数为80,垂直采样点数为300。
模型为一个三层水平层状均匀介质,其速度分别为1000m/s 、2000m/s 、3000m/s ,模型如下:
图1 模型示意图 图2 地震子波 3. 结果 2000m/s 1000m/s
3000m/s
图2 各层界面的地震响应。
子波整形提高合成地震记录质量单刚义,韩立国,张丽华,潘保芝吉林大学地球探测科学与技术学院,长春 130026 摘要:阐述合成地震记录的制作原理,利用测井数据计算出反射系数。
根据地震波传播机制,通过三维地震数据统计得到子波,由此计算出原始合成记录。
再根据过井剖面建立褶积模型,对子波进行整形,得到一个随时间变化的地震子波,最终得到与地震数据的波组特征相吻合的合成地震记录。
为精细储层描述打下了坚实的基础。
关键词:合成地震记录;褶积模型;子波整形;反射系数;深时转换;子波提取中图分类号:P 631.4 文献标识码:A收稿日期:2006207211作者简介:单刚义(19702),男,山东莱阳人,硕士研究生,主要从事地震资料的处理和解释研究工作,E 2m ail :ganyi 2shan@yahoo .com .cn .Usi ng W avelet Shap i ng to I m prove Syn thetic Se is m ogram Qua l itySHAN Gang 2yi ,HAN L i 2guo ,ZHAN G L i 2hua ,PAN B ao 2zh iColleg e of GeoE xp lora tion S cience and T echnology ,J ilin U n iversity ,Chang chun 130026,Ch inaAbstract :Syn thetic seis m ogram is the b ridge that connects seis m o logic data w ith w ell logging data ,and is also the basis of the in terp retati on of structu re and reservo ir’s seis m ic litho logy .T he accu racy of syn thetic seis m ogram directly affects the accu rate calib rati on of seis m ic geo logic ho rizon and the accu racy of reservo ir’s litho logy in terp retati on .T h rough th is p aper illu strates the p rinci p le of syn thesizingseis m ogram .T he reflecti on cofficien ts are ob tained w ith the w ell logging data .B y the p rinci p le of seis m ic w aves ,the w avelet is ex tracted from the seis m ic data of th ree di m en si on ,and the raw syn thetic seis m ogram is fin ished .T h rough bu ilding convo lu ti on m odel ,th is paper shapes w avelet ex tracted from seis m ic data ,th is w avelet varies w ith the ti m e .A t last ,the seis m ic litho logy of it is as the sam e of seis m ic data .T hu s it lays a strong foundati on fo r reservo ir fine descri p ti on .Key words :syn thetic seis m ogram ;convo lu ti on m odel ;w avelet shap ing ;reflecti on coefficien t ;dep th -ti m e conversi on ;ex tracti on of w avelet0 引 言测井和地震等资料的有机结合与综合运用是油藏描述中研究地质构造、岩性、岩相、油藏类型及预测有利含油气区等问题的基础。
合成地震记录业务流程一、准备工作。
咱得先把相关的数据都找齐喽。
比如说,测井数据那是相当重要的。
就像我们找宝藏得有个地图一样,测井数据就是我们合成地震记录的地图。
这里面包括声波测井曲线、密度测井曲线等。
这些数据就像是一个个小零件,缺了哪个都不行。
而且呀,我们还得确保这些数据的准确性,如果数据错了,那就好比做菜的时候盐当成了糖,做出来的东西肯定不对味。
另外呢,我们还需要有一些地质分层信息,这个就像是房子的框架结构,能让我们清楚地知道不同地层的情况,知道在哪个地层该怎么操作。
二、选择合适的子波。
子波就像是合成地震记录的画笔。
有好多不同类型的子波可以选呢。
我们得根据实际的地质情况和研究目的来挑。
如果是比较简单的地层结构,可能选个简单点的子波就够用啦。
但要是地层情况很复杂,就像一个超级复杂的迷宫一样,那我们就得找个功能强大、能适应复杂情况的子波。
这时候就得花点心思去对比不同的子波,看看哪个画出来的“画”(也就是合成的地震记录)最符合我们对这个地下情况的预期。
三、计算反射系数。
这一步就像是在做数学题,不过是很有趣的那种。
我们要根据前面准备好的测井数据,像声波和密度这些,来计算反射系数。
反射系数就像是镜子的反射率一样,它能告诉我们地震波在不同地层界面上反射的情况。
这个计算可不能马虎,要是算错了,那合成出来的地震记录就会像一个歪歪扭扭的积木塔,一点都不稳定也不准确。
我们得仔仔细细地按照公式来算,就像小心翼翼地搭积木一样,一块都不能搭错。
四、合成地震记录。
好啦,前面的工作都做好了,就到了最激动人心的合成地震记录这一步啦。
我们把选好的子波和计算好的反射系数放在一起,就像把颜料和画笔放在一起准备画画一样。
然后通过一些算法,让它们相互作用,就像魔法一样,一个地震记录就慢慢合成出来了。
这时候我们就像一个小魔法师,看着自己的作品一点点呈现出来。
不过呢,这时候还不能掉以轻心,我们还得检查这个合成出来的地震记录是不是合理。
五、验证与调整。
第一章时间序列分析基础一维傅里叶变换首先观察一个实验。
将弹簧的一端固定并悬垂,另一端挂一重物。
向下拉重物使弹簧拉伸某一距离,比如说0.8个单位,使其振动。
现假定弹簧是弹性的,那么它将无休止地上下运动。
若将运动起始的平衡位置定为时间零,那么重物的位移量将随着时间函数在极限[+0.8—-0.8]之间变化。
如果有一装置能给出位移振幅随时间函数变化的轨迹,就会得到一条正弦波曲线。
其相邻两峰值间的时间间隔为0.08秒(80毫秒)。
我们称它为弹簧的周期,它取决于所测弹簧刚度的弹性常数。
我们说弹簧在一个周期时间内完成了一次上下振动。
在1秒的观测时间内记下其周期数,我们发现是12.5周,这个数被称为弹簧振动的频率。
你一定会注意到,1/0.08=12.5,这就是说频率为周期的倒数。
我们取另一个刚性较大的弹簧,并重复上面的实验。
不过这次弹簧的振幅峰值位移为0.4个单位。
它的运动轨迹所显示的是另一条正弦曲线。
量其周期和频率分别为0.04秒和25周/秒,为了记下这些测量结果,我们做每个弹簧峰值振幅与频率的关系图,这便是振幅谱。
现在取两个相同的弹簧。
一个弹簧从0.8个单位的峰值振幅位移开始松开,并使其振动。
这时注意弹簧通过零时平衡位置的时间,就在它通过零时的一刹那,请你将另一弹簧从0.8个单位的同样峰值振幅位移处松开。
这样由于起始的最大振幅相同,所以两个正弦时间函数的振幅谱也应该一样。
但肯定两者之间是有差别的,特别是当第1个正弦波达到峰值振幅时,另一个的振幅为零。
两者的区别为:第2个弹簧的运动相对于第1个弹簧的运动有一个等于四分之一周期的时间延迟。
四分之一周期的时间延迟等于90°相位滞后。
所以除振幅谱之外,我们还可以作出相位延迟谱,至此,这个实验做完了。
那么我们学到了什么呢?这就是弹簧的弹性运动可以用正弦时间函数来描述,更重要的是,可以用正弦波的频率、峰值振幅及相位延迟来全面地描述正弦波运动。
这个实验告诉我们弹簧的振动是怎样随时间和频率函数变化的。
地震子波波形显示及一维地震记录合成一、实验目的1、认识地震子波(以雷克子波为例),对子波有直观的认识。
2、利用线性褶积公式合成一维地震记录。
二、实验内容1、雷克子波:零相位子波源:()()t f e t w m t f m πγπ2cos 22/2-= 程序:fm=30;r=3;t=0.002;for n=1:200w(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n);endfigure(1),plot(w);图像(1):图(1)最小相位子波:()()t f e t w m t f m πγπ2sin 22/2-= 程序:fm=10;r=3;t=0.002;for n=1:200w(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*sin(2*pi*fm*t*n);endfigure(1),plot(w);图像:图(2)其中m f 代表子波的中心频率,γ代表子波宽度,随着γ的增大,子波能量后移,当γ=7时,最小相位子波可视为混合相位子波。
因而若将最小相位子波中γ=7,则可以间接地得到混合相位子波的图形为:图(3)从图(3)中可以看出,经过改变 的值后,整个雷克子波的能量发生了明显的后移;从一定程度上可以反映最小能量子波、混合能量子波、最大能量子波的图像上的区别是能量集中区域上的不同。
2、利用线性褶积公式合成一维地震记录(1)利用零相位雷克子波作为震源合成一维地震记录:程序:fm=30;r=3;t=0.002;for n=1:200w(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n);endfor j=1:500r(j)=0;endr(100)=1.0;r(200)=-0.7;r(300)=0.5;r(400)=0.4;r(500)=0.6; 对地层界面的反射系数进行设定for j=1:699f(j)=0;for i=1:500if(j-i>0&&j-i<200)f(j)=f(j)+r(i)*w(j-i); 褶积的主要过程endendendfigure(2),plot(f);图像:图(4)在这个一维地震记录图像上可以看出,除了100、200、300、400、500这些存在反射系数的点上发生了变化,在这些点的周围也发生了些许的变化。
地震勘探原理实验一地震子波波形显示及一维地震合成记录姓名: 学号:专业:地球物理勘察技术与工程 级 一、实验目的1. 认识子波,对子波的波形有直观的认识。
(名词:零相位子波,混合相位子波,最小相位子波;了解子波的分辨率与频宽的关系;) 2. 利用褶积公式合成一维地震记录。
二、实验步骤 1. 雷克子波()()))(21(22t f et r m t f m ππ-=- 零相位子波())2sin()ln(222t f et w m n t f m π-= (最小相位子波) n= m1/m2为最大波峰m1和最大波谷m2之比()())2cos(log *22xw t f et w m mt f m +=-π 钟型子波 xw 为初相m 为时间域主波峰与次波峰之比w(t)=exp(-2*Fm^2*t^2*ln(n))*sin(T-2*pi*Fm*t) n=m1/m2 最大相位子波(最大相位子波请同学们自己查找相关文献完成,非必须完成)其中f代表子波的中心频率, t =i*dt,dt为时间采样间隔,i为时间m离散点序号; 这里可以为f = 10,25,40,100 Hz等,采样间隔dt=0.002m秒,i为0~256;2.根据公式编程实现不同频率的零相位子波的波形显示;不同中心频率的零相位子波图f = 25:mf = 100:m3.其地质模型为:设计反射系数)(n r (n=512),n 为地层深度,其中0.1)100(=r ,为第一层介质深度;7.0)200(-=r ,为第二层介质深度;5.0)300(=r ,为第三层介质深度;4.0)400(=r ,为第四层介质深度;6.0)450(=r ,为第五层介质深度;其它为0。
地震波在介质中传播,当到达介质分界面时,发生反射和透射,反射波被检波器接受,生成地震记录。
反射系数表示地震波在两层介质分界面的能量重新分配,如r(100)=1.0,表示地震波入射到分界面时,只有一种波,反射纵波(或反射横波)。
反射系数不为1.0时,表示当地震波入射到分界面时,产生两种反射波。
反射系数为正,表示反射波相位与入射波相位相差2π;反射系数为负,表示反射波相位与入射波相位相差π,存在半波损失。
4. 应用褶积公式∑=-=*=Nm m n w m r n w n r n f 1)()()()()(合成一维地震记录,并图形显示;应用褶积公式求f (n )的程序为:#include<stdio.h> #include<math.h> #define PI 3.1415926 #define FM 100 void main() {double fac(double x[],double y[],double z[],int m,int n); FILE *fp;int i,j,x;double W,dt=0.002,t,a[256];double b[512]={0};double r[512]={0};r[100]=1.0;r[200]=-0.7;r[300]=0.5;r[400]=0.4;r[450]=0.6;fp=fopen("Date.txt","w+");printf("please input x:\n");scanf("%d",&x);for(i=0;i<256;i++){t=i*dt;if(x==1)W=exp(-2*FM*FM*t*t*log(2))*sin(2*PI*FM*t); else if (x==2)W=(1-2*pow(PI*FM*t,2))*exp(-pow(PI*FM*t,2)); else if (x==3)W=exp(-FM*FM*t*t*log(2))*cos(2*PI*FM*t+PI/4);a[i]=W;}fac(r,a,b,512,256);for(j=0;j<512;k++){fprintf(fp,"%f\n",b[j]);}}double fac(double x[],double y[],double z[],int m,int n){int i,j;for(i=0;i<=m+n-1;i++){double sum=0.0;for(j=0;j<=m;j++){if(i-j>0&&i-j<=256)sum+=x[j]*y[i-j];}z[i]=sum;}}三、实验结果一维反射系数序列图形显示为:零相位子波与反射系数褶积后的地震记录图形显示:f = 25:mf = 100:m最小相位子波与反射系数褶积后的地震记录图形显示:f = 25:mf = 100:m混合相位子波与反射系数褶积后的地震记录图形显示:f = 25:mf = 100:m最大相位子波与反射系数褶积后的地震记录图形显示:f = 25:mf = 100:m零相位振幅图形显示:f = 25:mf = 100:m零相位幅角图形显示:f = 25:mf = 100:m最小相位振幅图形显示:f = 25:mf = 100:m最小相位幅角图形显示:f = 25:mf = 100:m混合相位幅角图形显示:f = 25:mf = 100:m混合相位振幅图形显示:f = 25:mf = 100:m最大相位幅角图形显示:f = 25:mf = 100:m最大相位振幅图形显示:f = 25:mf = 100:m四、实验分析根据所学知识对实验结果进行分析;地震子波由震源激发,在地层中传播,因为在沉积地层中,每层介质的物理性质不相同,从而使得地震波的传播速度也不相同。
当地震波传播到两层介质的分界面时,会发生反射,由于每层介质的反射系数不同,所以反射波的能量也不相同,检波器接收到不同时刻的、不同能量的反射波,形成一个地震记录。
由合成地震记录中可以看出,最小相位子波相对零相位子波来说是相位滞后的,能量延迟的,但两者为同一家族的子波。
合成地震记录中横坐标为时间,纵坐标为振幅。
每一时刻的值由m个值的和组成,m为反射系数r(n)的长度,整个地震记录由m+n-1个时刻的值组成。
对于零相位的地震记录来说,当r(m)=1.0时,即j=100时,i=100时,w(i-j)=1.0,是能量最大的,即w(0)=1.0。
同理,当n=200,、300、400、450时,能量也是最大的。
对于最小相位的地震记录来说,当r(m)=1.0时,即j=100时,但i=100时,w(i-j)不是最大能量的,即最大能量不是在w(0)出现,而是延迟出现。
同理,当n=200、300、400、450时,能量也不是最大的,而是要延迟出现。
由振幅图及幅角图可知,零相位子波能量聚集在首部,开始时就具有最大能量,无积累过程,当振幅最大时,相位为零,即此时波的振幅为实数,达到最大值;最小相位子波能量聚集在序列首部,是最小能量延迟的,信号随时间增大而减小,当振幅最大时,相位不为零,是非零相位的,相对零相位子波来说,最大能量是延迟的;混合相位子波的能量聚集在序列中部,是混合能量延迟的;最大相位子波能量聚集在后部。
最大相位子波和混合相位子波的信号信号不随时间增大而减小。
五、附:源程序代码#include<stdio.h>#include"13KFFT.C"#include<math.h>#define PI 3.1415926#define FM 100void main(){double fac(double x[],double y[],double z[],int m,int n);FILE *fp,*fpr,*fpre,*fpi,*fpamp,*fpha;int i,j,x;doubleW,dt=0.002,t,a[256],pr[512],pi[512]={0.0},fr[512],fi[512],amp[512],p ha[512];double b[512]={0};double r[512]={0};r[100]=1.0;r[200]=-0.7;r[300]=0.5;r[400]=0.4;r[450]=0.6;fp=fopen("褶积结果.txt","w+");fpr=fopen("反射系数.csv","w+");fpre=fopen("实部.txt","w+");fpi=fopen("虚部.txt","w+");fpamp=fopen("振幅.csv","w+");fpha=fopen("相位.csv","w+");for(i=0;i<512;i++){fprintf(fpr,"%f\n",r[i]);}fclose(fpr);printf("please input x:\n");scanf("%d",&x);for(i=0;i<256;i++){t=i*dt;if(x==1)W=(1-2*pow(PI*FM*t,2))*exp(-pow(PI*FM*t,2));else if (x==2)W=exp(-2*FM*FM*t*t*log(2))*sin(2*PI*FM*t);else if (x==3)W=exp(-2*pow(FM*t,2)*log(2))*sin(0.512-2*PI*FM*t);else if(x==4)W=exp(-FM*FM*t*t*log(2))*cos(2*PI*FM*t+PI/4); a[i]=W;}fac(r,a,b,512,256);for(j=0;j<512;j++){fprintf(fp,"%f\n",b[j]);}for(i=0;i<512;i++){pr[i]=b[j];}for(i=0;i<512;i++){pr[i]=fr[i+127];}kfft(pr,pi,512,9,fr,fi,0,1);for(i=0;i<512;i++){fprintf(fpre,"%e\n",fr[i]);fprintf(fpi,"%e\n",fi[i]); fprintf(fpamp,"%e\n",pr[i]);fprintf(fpha,"%f\n",pi[i]);}fclose(fpre);fclose(fpi);fclose(fpamp);fclose(fpha);}double fac(double x[],double y[],double z[],int m,int n) {int i,j;for(i=0;i<=m+n-1;i++){double sum=0.0;for(j=0;j<=m;j++){if(i-j>0&&i-j<=256)sum+=x[j]*y[i-j];}z[i]=sum;}}。