基于MATLAB的地震正演模型实现.docx

1、基于MATLAB的地震正演模型实现基于MATLAB的地震正演模型实现贾跃玮;杨锐【摘 要】人工合成地震正演模型是进行三维模型计算的基础,在地震勘探领域具有重要意义.针对地震勘探原理,作者运用MATLAB强大数学计算和图像可视化功能,对一个三层介质模型制作二维人工合成地震记录.文章说明地震记录形成的物理机制,介绍了地质模型的构造及参数选择,针对具体地质模型制作合成地震记录,成功验证了褶积模型原理.【期刊名称】计算机与数字工程【年(卷),期】2009(037)007【总页数】4页(P132-135)【关键词】地震;MATLAB;正演【作 者】贾跃玮;杨锐【作者单位】中国地质大学地下信息探测技术与仪

2、器教育部重点实验室,北京,100083;川庆钻探工程有限公司地质勘探开发研究院,成都,610059【正文语种】中 文【中图分类】工业技术【文献来源】总第 237 期2009 年第 7 期 计算机与数字工程 Computer&DigitalEngineering Vol.37No.7132基 于 1/IATLAB的 地 震 正 演 模 型 实 现+贾跃玮 杨 孙 （中国地质大学地下信息探测技术与仪器教育部重点实验室” 北京100083 ） （川庆钻探工程有限公司地质勘探开发研究院” 成都 610059 ）摘 要 人工合成地震正演模型是进行三维模型计算的基础 ，在地震勘探领域具有重要意 义 。 针

3、对地震勘探原理 ，作者运用 MATLAB强大数学计算和图像可视化功能 ，对一个三层介质模型制作二维人工合成地震记录 。 文章说明地震记录形成的物理机制 ，介绍 了地质模型的构造及参数选择 ，针对具体地质模型制作合成地震记录 ，成功验证 了褶积模 型原理 。 关键词 地震 MATI 。AB正演中图分类号P631.4ASeismicForward ModelingBasedonMATLAB Jia Yuewei YangRuiz) (Geo-detectionLab,Ministry of EducationChinaUniversityof Geosciences.Beijing 100083)

4、 (GeologicalExplorationandDevelopmentResearchInstituteof SichuanPetroleumAdministration2) ,Chengdu 610059) Abstract A syntheticseismicforwardmodelis thebaseof 3Dmodelingcalculation,andhasprofoundmeaningin seismicexploration.Basedontheprinciplesofseismic exploration,theauthorappliedgreatmathematiccal

5、culationandfigural visualization of MATLABtobuild asyntheticseismic recordingof athree-layermodel.Thephysicalprincipleofseismic explorationisintroducedatthe firstpart of the article.Then,thetectonic structureof the geophysicalmodelandparametersaredelineated.In thelast partof thepaper,the methodof sy

6、ntheticseismic recordsbasedonthe specifiedge-ophysicalmodelispresented. Keywordsseismic,MATLAB,forward modeling ClassNumber P631.4 1 引言r地震勘探是利用地下介质弹性和密度的差异 ， 通过观测和分析大地对人工激发地震波 的响应 ，推 断地下岩层 的性质和形 态 的地球物理方法 。 地震 勘探是钻探前勘测石油与天然气资源 的重要手段 ， 在煤 田和工程地质勘查 、 区域地质研究和地壳研究 等方面 ，也得到广泛应用 。人工合成二维地震模型 记 录是各种 复杂地震 模型

7、正演计算的基础 ，是对地震勘探经典理论的忠 实实现 。 在实际工作中，针对具体地质构造进行二 维地震模拟能够有效帮助地球物理工作者在地震剖面上 识 别各 种 地 质 现象 。 MATLAB 是 Math- works公司于 1984 年推 出 的一 款数学 软件 ，是用 于工程科学计算 的高效率编 程语 言 23 。MAT- LAB不仅在数值科学计算上 独 占鳌头 ， 而且 图形 可视化和 图像处理能力也毫不逊色n6 ，特别适合人工 合 成 地 震 记 录 的 快 速 实 现 。 因 此 ， 本 文 在 MATI 。AB 环境下设计 了一个 三层地质模 型 ，并对该模型模拟了地震记录 ， 旨在

8、可视化地观察地震波场记录特征并验证地震褶积模型 。 2地震记 录形成 的物 理机制在地震记录上看到 的波形 是地震子波叠加 的 *收稿 日期 ：2009 年 3 月 23 日，修回 日期 ：2009 年 4 月 21 日作者简介 ：贾跃玮 ，男 ，硕士研究生 ，研究方向：地震资料解释。 杨锐 ，研究方向：勘察地球物理。总第 237 期 2009 年第 7 期计算机与数字工程 Computer&DigitalEngineering Vol.37No.7 132贾跃玮杨孙（”北京 100083 ）成都610059）摘要人工合成地震正演模型是进行三维模型计算的基础 ，在地震勘探领域具有重要意 义 。

9、 针对地震勘探原理 ，原理 。关键词地震P631.4 SeismicForward ModelingBasedonMATLAB Yang Ruiz) (Geo-detectionLab,Ministry of EducationChinaUniversityof Geosciences.Beijing (GeologicalExplorationandDevelopmentResearchInstituteof SichuanPetroleumAdministration2) ,Chengdu syntheticseismicforwardmodelis thebaseof 3Dmodelin

10、gcalculation,andhasprofoundmeaningin seismicexploration.Basedontheprinciplesofseismic exploration,theauthorappliedgreatmathematiccalculationand figural visualization of MATLABtobuild asyntheticseismic recordingof athree-layermodel.Thephysicalprincipleof seismic explorationisintroducedatthe firstpart

11、 of the article.Then,thetectonic structureof the geophysicalmodeland parametersaredelineated.In thelast partof thepaper,the methodof syntheticseismic recordsbasedonthe specifiedge- ophysicalmodelispresented. Key words seismic,MATLAB, forward ClassNumber 1引言 r地震勘探是利用地下介质弹性和密度的差异 ，通过观测和分析大地对人工激发地震波 的响

12、应 ，推断地下岩层 的性质和形 态 的地球物理方法 。 地震勘探是钻探前勘测石油与天然气资源 的重要手段 ，在煤 田和工程地质勘查 、 区域地质研究和地壳研究等方面 ，也得到广泛应用 。人工合成二维地震模型 记 录是各种 复杂地震模型正演计算的基础 ，是对地震勘探经典理论的忠实实现 。 在实际工作中，针对具体地质构造进行二维地震模拟能够有效帮助地球物理工作者在地震公司于 1984 年推 出 的一 款数学 软件 ，是用于工程科学计算 的高效率编 程语 言 23。不仅在数值科学计算上 独 占鳌头 ， 而且 图形可视化和 图像处理能力也毫不逊色n6，特别适合 MATI。AB环境下设计 了一个 三层地

13、质模 型 ，并对第 37 卷(2009) 第 7 期结果 ，从地下许多反射界面发生反射时形成的地震 子波 ，振幅大小 决定 于 反 射界 面 反 射 系数 的绝 对 值 ，极性的正负决定 于反 射 系数 的正 负 ， 到达时 间 的先后取决于界面深度 和覆盖层 的波速 1 。若地震子波波形用 S(t) 表示 ，反 射 系数是 双程垂直反 射旅行时 f 的函数 ，用 R(t) 表示 ，地震 记录 厂 (z)形 成的物理过程在数学上就可表示为 ： (f) 一 S(f)*R(f)一S(r)R(t-r)dr地震子波和反射系数资料常常不易取得 ， 因此计算时常做这样一些假设 ： 1) 地质模型 的建立是

14、来 自大量观察实 际地质结构的经验性归纳总结 。 2) 为了模型建立和计算过程 中突出理论数值 ， 去除 了一些干扰 因素 ，对一切 衰减 、 噪 声都不进行 考虑 。 3) 地层在横 向上均匀 ，纵向上是 由大量具有不同弹性性质的薄层构成 。 4) 地震子波以平面波形式垂直入射到界面 ，各 薄层 的反射子波与地震子波形状相 同， 只是振幅及 极性不 同。 5) 所有波的转换 、 吸收及绕射等能量损失都不考虑 。基于以上假设条件进行地震记 录合成就必须 已 知地震子波 以及地层 的反射系数 ，而反射系数又 主要 由地层 的波阻抗反映 ，所以必须首先获取地层 的速度和密度资料 。 速度资料可通过

15、 连续 速度测 井获得 ，密度资料可从密度测 井获得 ，得不 到密度 资料时 ，可近似假定密度不变 ， 以速度 曲线代替波 阻抗 曲线来计算反射系数 。 加德纳根据实际 资料 提 出了一个 由速度推算密度 的经验公式 ：10=0.23V25( 速度单位 ：ft/s)或 P=0.31V25( 速度单位 ：m/s)没有速度测井资料时 ，若有 电阻率测井 曲线 ，则可用法斯特公式 ： V =KH告 R: 其中，V是速度 ，K 是一个与岩石性质有关 的参数 ， R是 电阻率 ，H 是深度 。已知地震子波是合成地震记 录 的一个很重要 的前提条件 。 在 已 有 的地震记录上选取地震子波 的具体方法如下

16、 ： 1)在地震记录上识别 出单波 ，做 出单波波形 ， 再反复试验 ，检查找 出符合实际的子波 。 2) 根据 已总结 出的地震子波特点 ，采取具有特 殊数学表达式的波形表示 ，如雷克子波等。3) 使用非炸药震源时记录的震源子波波形 。4) 利用实际地震记录 ，采取数字处理方法在一 定的假设条件下求取地震子波 。 5) 在拥有井中观测初至纪录时 ，可考虑使用初至波做子波波形 。 6) 在拥有声波测井资料和井旁地震 记录 X(t) 时 ，反射系数曲线 R(t) ，地震子波 S(t) ，可 由 X( 厂 )=S( )*R(f) 得子波 的谱 S(f)=X( 厂 )/R( ) ， 再对上式作反傅氏

17、变换得地震子波波形 S(t) 。 3地质模型 的建 立地质及地球物理学研究表明 ，地表之下 的地质 结构是极其复杂 的。 地质构造是指地壳 中的地层 在岩层地壳运 动 的作用 下发生变形 与变位而遗 留 下来的形态本 文 中采用 的是一 个 三层水平层状均匀介质模型 ， 即似没有发育在三个不 同时期 的 水平层状沉积岩层 。 同时 ，为 了体现地震勘探在石 油 、 天然气勘探领域 的应用效果 ，作者结合实 际地 质资料模拟 了地下储层含有石油的情况 。 一般情况下 ，随着埋藏深度越来越大 ，岩石的地 震波传播速度会越来越大。 深层的岩石地震波传播 速度通常要大于浅层岩石的传播速度 。 而流体的

18、传 播速度较岩石的传播速度要小很多。 比如，砂岩的 地震波传播速度一般在 1800m/s4000m/s，石油 的传播速度一般为 1300m/s 1400m/s。 在模型建 立的过程中，作者参照 了实际地质构造及各种岩石 的速度资料 ，确定 了三层水平层状均匀介质的速度 及深度参数。 储层含油地质模型如图 1 所示 ： 第一层 H-IOOOm V=1500m/s 第二层 H=1200m含油层 H=1200m V=2000m/s V=1300m/s V=2000m/s 第三层 H=1400m V=2500m/s 图 1 含油地质模型 第 一 层 H=1000m V=1500m/s V=2000m/

19、s V=2500m/s 在该地 质模 型 中，作 者设 计了一个含油层 ，并使该含油层 包裹在砂岩环境中，各层的速度 及深度参数如模型 中所示 。 此 外 ，为了更加直观地显示地震勘 探的效果 ，本文提供了一个参照 图2对比地质模型 模型 ，该模型中没有设计含油速第37卷(2009) 第 7 期结果 ，从地下许多反射界面发生反射时形成的地震子波 ，振幅大小 决定 于 反 射界 面 反 射 系数 的绝 对值 ，极性的正负决定 于反 射 系数 的正 负 ， 到达时 间的先后取决于界面深度 和覆盖层 的波速 1震子波波形用 S(t) 表示 ，反 射 系数是 双程垂直反射旅行时 f 的函数 ，用 R(

20、t) 表示 ，地震 记录 厂 (z)形成的物理过程在数学上就可表示为 ：(f) 一 S(f)*R(f)一 S(r)R(t-r)dr 2) 为了模型建立和计算过程 中突出理论数值 ，去除 了一些干扰 因素 ，对一切 衰减 、 噪 声都不进行 4) 地震子波以平面波形式垂直入射到界面 ，各薄层 的反射子波与地震子波形状相 同， 只是振幅及极性不 同。基于以上假设条件进行地震记 录合成就必须已 知地震子波 以及地层 的反射系数 ，而反射系数又主要 由地层 的波阻抗反映 ，所以必须首先获取地层的速度和密度资料 。 速度资料可通过 连续 速度测井获得 ，密度资料可从密度测 井获得 ，得不 到密度资料时

21、，可近似假定密度不变 ， 以速度 曲线代替波阻抗 曲线来计算反射系数 。 加德纳根据实际 资料提 出了一个 由速度推算密度 的经验公式 ： 10=0.23V25( 速度单位 ：ft/s)告R:其中，V已知地震子波是合成地震记 录 的一个很重要的前提条件 。 在 已 有 的地震记录上选取地震子波的具体方法如下 ：在地震记录上识别 出单波 ，做 出单波波形 ，再反复试验 ，检查找 出符合实际的子波 。 2) 根据 已总结 出的地震子波特点 ，采取具有特殊数学表达式的波形表示 ，如雷克子波等。 3) 使用非炸药震源时记录的震源子波波形 。 4) 利用实际地震记录 ，采取数字处理方法在一定的假设条件下

22、求取地震子波 。 6) 在拥有声波测井资料和井旁地震 记录 X(t)时 ，反射系数曲线 R(t) ，地震子波 S(t) ，可 由 X( 厂 ) =S()*R(f)得子波 的谱 S(f)=X( 厂 )/R( )，再对上式作反傅氏变换得地震子波波形 S(t) 。地质及地球物理学研究表明 ，地表之下 的地质结构是极其复杂 的。 地质构造是指地壳 中的地层在岩层地壳运 动 的作用 下发生变形 与变位而遗 留下来的形态状均匀介质模型 ， 即似没有发育在三个不 同时期 的水平层状沉积岩层 。 同时 ，为 了体现地震勘探在石油 、 天然气勘探领域 的应用效果 ，作者结合实 际地质资料模拟 了地下储层含有石油

23、的情况 。一般情况下 ，随着埋藏深度越来越大 ，岩石的地震波传播速度会越来越大。 深层的岩石地震波传播速度通常要大于浅层岩石的传播速度 。 而流体的传播速度较岩石的传播速度要小很多。 比如，砂岩的地震波传播速度一般在 1800m/s4000m/s，石油的传播速度一般为 1300m/s 1400m/s。 在模型建立的过程中，作者参照 了实际地质构造及各种岩石的速度资料 ，确定 了三层水平层状均匀介质的速度及深度参数。 储层含油地质模型如图 1 所示 ： V=1500 m/s V=2000 V=1300 V=2500图1含油地质模型 V=1 500在该地 质模 型 中，作 者设计了一个含油层 ，并

24、使该含油层包裹在砂岩环境中，各层的速度及深度参数如模型 中所示 。 此外 ，为了更加直观地显示地震勘探的效果 ，本文提供了一个参照134贾跃玮等 ：基于 MATLA13 的地震正演模型实现第 37 卷度突变层 ，在地震模拟记录参数不边 的条件下用来 比照含油模型的效果 。 模型构造及参数如图 2 。 4地震模型正演根据褶 积理论 ，结合地质模 型 ，作 者在 MAT- I 。AB 环境中编写程序实现 了二维地震正演。 首先 ，为了最大可能的与实际情况相符合，地震模型 中使用的子波是稳定可实现的子波 8，如图 3 。图 3 子波 该子波是最小相位子波 ，有时称为前载子波 ，其能量集中在整个波形

25、的前端。 由于大多数脉 冲 地震震源 （ 如炸药震源 ） 产生 的原始脉 冲是接近最 小相位的，因此在地震正演模型 中的地震子波选取 一般都选择最小相位型 子波 7二维地震正演模型实现的主要程序及关键步骤注释如下 ： n=5000 ； 采样点数 dx=50； 道间距dt 一0.002； 采样间隔nl=3； 界面数m- 80； 道数v0= 1000150020002500 ；第 1 层速度vl= 1000150013002500 ；第 2 层速度v2= 1000150020002500 ；第 3 层速度h(l，1)=800； 7 h(2，1)=1200： h(3，1)=1800该模型为各层平行

26、，无倾角 xmax=(m/2+1)*dx; ymax=n*dt； 坐标范围p=2;图像中显示数字以下部分为子波采样a0=200;f=20 ；nw-60 b=30; 子波参数 tt=O:dt:(nw-l)* dt; wb=a0*sin(2pi*f*tt).*exp( -b*tt) ； 视速度hvaOh(3,:).*vO(1:n1); hvbO=h(3,:). /vO(l:nl); hva1一h(3.:).*v1(1:nl);hvbl=h(3,:)./vl(l:nl);hva2一 h(3,:).*v2(1:n1);hvb2=h(3,:)./v2(l:nl); for i=2:nl hva0(i)一

27、hva0(i)+hva0(i-1):hvbO(i)=hvbO(i)+hvbO(i-l);hval(i)=hval(i)+hval(i-l);hvbl(i)=hvbl(i)+hvbl(i-l):hva2(i) =hva2(i)+hva2(i-l):hvb2(1)=hvb2(i)+hvb2(1-1);end hva0 sqrt(hva0./hvbo);hval=sqrt(hval./hvbl);hva2一 sqrt(hva2./hvb2);hva0(1) 一v0(1)：hval(l)=vl(l) ： hva2(1) 一v2(1) ；以上程序段使用循环控制各个层速度变化 for i=l:nl设定反射

28、面深度 forj=2:mh(i,j)=h(i,j-l)end endmk=round(m/2) ；炮检距 hvbO(mk)=0;hvbl(mk)=0;hvb2( mk)=o ：hvbl=(l:m)*dx;hvbl= (l:m)*dx;hvb2= (l:m)*dx;for i=l:nl反射系数 r0(i) 一 3*(v0( i+l)-v0(i) (v0(i+l)+v0(i); r1(i) 一3*(v1(i+1)一 v1(i) (vl(i+1)+vl(i); r2(i) 一3*(v2(i+1)一v2(i) (v2(i+1)+v2(i); endx=zeros(n,m) b=pi/180； 划分网隔

29、for i=1 ：nl一次反射波 forj=1:80ifj30& 。j 一50 ；z=2* h(i,j)/hval(i) ； t round(z/dt) ：x(t，j)一x(t ，j)+r1(i) ；else z=2*h(i,j)/hva2(i);t= round(z/dt) ：x(t，j)一x(t ，j)+r2(i) ；end卷度突变层 ，在地震模拟记录参数不边 的条件下用来比照含油模型的效果 。 模型构造及参数如图 2 。 I环境中编写程序实现 了二维地震正演。 首先 ，3子波该子波是最小相位子波 ，有时称为前载子波 ，其能量集中在整个波形 的前端。 由于大多数脉 冲地震震源 （ 如炸药震

30、源 ） 产生 的原始脉 冲是接近最小相位的，因此在地震正演模型 中的地震子波选取一般都选择最小相位型 子波 7 n= 5000 ；采样点数 dx= 50；道间距 dt 一0.002采样间隔 nl=3界面数 m- 80；道数 v0= 1000150020002500 ；层速度 vl= 1000150013002500 ； v2= 1000150020002500 ； h(l，1)=8007 xmax= (m/2+1)* dx;坐标范围 p=2;以下部分为子波采样 a0=200;f=20nw- 60b=30;子波参数 wb= a0* sin(2pi*f*tt).*exp( -b*tt)视速度 hv

31、aOh(3,:).*vO(1:n1); hva1一h(3.:).*v1(1:nl); hvbl=h(3,:)./vl(l:nl); hva2一 h(3,:).*v2(1:n1); hvb2=h(3,:)./v2(l:nl); hva0(i)一hva0(i)+hva0(i-1): hvbO(i)=hvbO(i)+hvbO(i-l); hval(i)=hval(i)+hval(i-l); hvbl(i)=hvbl(i)+hvbl(i-l): hva2(i) =hva2(i)+hva2(i-l): hvb2(1)=hvb2(i)+hvb2(1-1); end hva0 sqrt(hva0./hvbo); hval=sqrt(hval./hvbl); hva2一 sqrt(hva2./hvb2); hva0(1) 一v0(1)： hval(l)=vl(l) forj=2:m h(i,j)=h(