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

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

《基于MATLAB的地震正演模型实现.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的地震正演模型实现.docx（7页珍藏版）》请在冰豆网上搜索。

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

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

贾跃玮;杨锐

【摘要】人工合成地震正演模型是进行三维模型计算的基础,在地震勘探领域具有重要意义.针对地震勘探原理,作者运用MATLAB强大数学计算和图像可视化功能,对一个三层介质模型制作二维人工合成地震记录.文章说明地震记录形成的物理机制,介绍了地质模型的构造及参数选择,针对具体地质模型制作合成地震记录,成功验证了褶积模型原理.

【期刊名称】《计算机与数字工程》

【年（卷）,期】2009（037）007

【总页数】4页（P132-135）

【关键词】地震;MATLAB;正演

【作者】贾跃玮;杨锐

【作者单位】中国地质大学地下信息探测技术与仪器教育部重点实验室,北京,100083;川庆钻探工程有限公司地质勘探开发研究院,成都,610059

【正文语种】中文

【中图分类】工业技术

【文献来源】

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

针对地震勘探原理，作者运用MATLAB强大数学计算和图像可视化功能，对一个三层介质模型制作二维人工合成地震记录。

文章说明地震记录形成的物理机制，介绍了地质模型的构造及参数选择，针对具体地质模型制作合成地震记录，成功验证了褶积模型原理。

关键词地震MATI。

AB正演中图分类号P631.4ASeismicForwardModelingBasedonMATLABJiaYuewei"YangRuiz）（Geo-detectionLab,MinistryofEducationChinaUniversityofGeosciences".Beijing100083）（GeologicalExplorationandDevelopmentResearchInstituteofSichuanPetroleumAdministration2）,Chengdu610059）AbstractAsyntheticseismicforwardmodelisthebaseof3Dmodelingcalculation,andhasprofoundmeaninginseismicexploration.Basedontheprinciplesofseismicexploration,theauthorappliedgreatmathematiccalculationandfiguralvisualizationofMATLABtobuildasyntheticseismicrecordingofathree-layermodel.Thephysicalprincipleofseismicexplorationisintroducedatthefirstpartofthearticle.Then,thetectonicstructureofthegeophysicalmodelandparametersaredelineated.Inthelastpartofthepaper,themethodofsyntheticseismicrecordsbasedonthespecifiedge-ophysicalmodelispresented.Keywordsseismic,MATLAB,forwardmodelingClassNumberP631.4／／1引言r地震勘探是利用地下介质弹性和密度的差异，通过观测和分析大地对人工激发地震波的响应，推断地下岩层的性质和形态的地球物理方法。

地震勘探是钻探前勘测石油与天然气资源的重要手段，在煤田和工程地质勘查、区域地质研究和地壳研究等方面，也得到广泛应用。

人工合成二维地震模型记录是各种复杂地震模型正演计算的基础，是对地震勘探经典理论的忠实实现。

在实际工作中，针对具体地质构造进行二维地震模拟能够有效帮助地球物理工作者在地震剖面上识别各种地质现象。

MATLAB是Math-works公司于1984年推出的一款数学软件，是用于工程科学计算的高效率编程语言‘2~3]。

MAT-LAB不仅在数值科学计算上独占鳌头，而且图形可视化和图像处理能力也毫不逊色n~6]，特别适合人工合成地震记录的快速实现。

因此，本文在MATI。

AB环境下设计了一个三层地质模型，并对该模型模拟了地震记录，旨在可视化地观察地震波场记录特征并验证地震褶积模型。

2地震记录形成的物理机制在地震记录上看到的波形是地震子波叠加的*收稿日期：

2009年3月23日，修回日期：

2009年4月21日作者简介：

贾跃玮，男，硕士研究生，研究方向：

地震资料解释。

杨锐，研究方向：

勘察地球物理。

总第237期2009年第7期计算机与数字工程Computer&DigitalEngineeringVol.37No.7132贾跃玮¨杨~孙（”北京100083）成都610059）摘要人工合成地震正演模型是进行三维模型计算的基础，在地震勘探领域具有重要意义。

针对地震勘探原理，原理。

关键词地震P631.4SeismicForwardModelingBasedonMATLABYangRuiz）（Geo-detectionLab,MinistryofEducationChinaUniversityofGeosciences".Beijing（GeologicalExplorationandDevelopmentResearchInstituteofSichuanPetroleumAdministration2）,Chengdusyntheticseismicforwardmodelisthebaseof3Dmodelingcalculation,andhasprofoundmeaninginseismicexploration.Basedontheprinciplesofseismicexploration,theauthorappliedgreatmathematiccalculationandfiguralvisualizationofMATLABtobuildasyntheticseismicrecordingofathree-layermodel.Thephysicalprincipleofseismicexplorationisintroducedatthefirstpartofthearticle.Then,thetectonicstructureofthegeophysicalmodelandparametersaredelineated.Inthelastpartofthepaper,themethodofsyntheticseismicrecordsbasedonthespecifiedge-ophysicalmodelispresented.Keywordsseismic,MATLAB,forwardClassNumber1引言r地震勘探是利用地下介质弹性和密度的差异，通过观测和分析大地对人工激发地震波的响应，推断地下岩层的性质和形态的地球物理方法。

地震勘探是钻探前勘测石油与天然气资源的重要手段，在煤田和工程地质勘查、区域地质研究和地壳研究等方面，也得到广泛应用。

人工合成二维地震模型记录是各种复杂地震模型正演计算的基础，是对地震勘探经典理论的忠实实现。

在实际工作中，针对具体地质构造进行二维地震模拟能够有效帮助地球物理工作者在地震公司于1984年推出的一款数学软件，是用于工程科学计算的高效率编程语言‘2~3]。

不仅在数值科学计算上独占鳌头，而且图形可视化和图像处理能力也毫不逊色n~6]，特别适合MATI。

AB环境下设计了一个三层地质模型，并对第37卷（2009）第7期结果，从地下许多反射界面发生反射时形成的地震子波，振幅大小决定于反射界面反射系数的绝对值，极性的正负决定于反射系数的正负，到达时间的先后取决于界面深度和覆盖层的波速‘1]。

若地震子波波形用S（t）表示，反射系数是双程垂直反射旅行时f的函数，用R（t）表示，地震记录厂（z）形成的物理过程在数学上就可表示为：

／’（f）一S（f）*R（f）一S（r）R（t-r）dr地震子波和反射系数资料常常不易取得，因此计算时常做这样一些假设：

1）地质模型的建立是来自大量观察实际地质结构的经验性归纳总结。

2）为了模型建立和计算过程中突出理论数值，去除了一些干扰因素，对一切衰减、噪声都不进行考虑。

3）地层在横向上均匀，纵向上是由大量具有不同弹性性质的薄层构成。

4）地震子波以平面波形式垂直入射到界面，各薄层的反射子波与地震子波形状相同，只是振幅及极性不同。

5）所有波的转换、吸收及绕射等能量损失都不考虑。

基于以上假设条件进行地震记录合成就必须已知地震子波以及地层的反射系数，而反射系数又主要由地层的波阻抗反映，所以必须首先获取地层的速度和密度资料。

速度资料可通过连续速度测井获得，密度资料可从密度测井获得，得不到密度资料时，可近似假定密度不变，以速度曲线代替波阻抗曲线来计算反射系数。

加德纳根据实际资料提出了一个由速度推算密度的经验公式：

10=0.23V'25（速度单位：

ft/s）或P=0.31V25（速度单位：

m/s）没有速度测井资料时，若有电阻率测井曲线，则可用法斯特公式：

V=KH告R~:

其中，V是速度，K是一个与岩石性质有关的参数，R是电阻率，H是深度。

已知地震子波是合成地震记录的一个很重要的前提条件。

在已有的地震记录上选取地震子波的具体方法如下：

1）在地震记录上识别出单波，做出单波波形，再反复试验，检查找出符合实际的子波。

2）根据已总结出的地震子波特点，采取具有特殊数学表达式的波形表示，如雷克子波等。

3）使用非炸药震源时记录的震源子波波形。

4）利用实际地震记录，采取数字处理方法在一定的假设条件下求取地震子波。

5）在拥有井中观测初至纪录时，可考虑使用初至波做子波波形。

6）在拥有声波测井资料和井旁地震记录X（t）时，反射系数曲线R（t），地震子波S（t），可由X（厂）=S（／’）*R（f）得子波的谱S（f）=X（厂）/R（／’），再对上式作反傅氏变换得地震子波波形S（t）。

3地质模型的建立地质及地球物理学研究表明，地表之下的地质结构是极其复杂的。

地质构造是指地壳中的地层在岩层地壳运动的作用下发生变形与变位而遗留下来的形态…本文中采用的是一个三层水平层状均匀介质模型，即似没有发育在三个不同时期的水平层状沉积岩层。

同时，为了体现地震勘探在石油、天然气勘探领域的应用效果，作者结合实际地质资料模拟了地下储层含有石油的情况。

一般情况下，随着埋藏深度越来越大，岩石的地震波传播速度会越来越大。

深层的岩石地震波传播速度通常要大于浅层岩石的传播速度。

而流体的传播速度较岩石的传播速度要小很多。

比如，砂岩的地震波传播速度一般在1800m/s～4000m/s，石油的传播速度一般为1300m/s～1400m/s。

在模型建立的过程中，作者参照了实际地质构造及各种岩石的速度资料，确定了三层水平层状均匀介质的速度及深度参数。

储层含油地质模型如图1所示：

┏━━━━━━━━━━━━━━━━━━━━━━━━━━┓┃第一层H-IOOOm┃┃V=1500m/s┃┣━━━━━━━━┳━━━━━━━━┳━━━━━━━━┫┃第二层H=1200m┃含油层H=1200mV=2000m/sV=1300m/sV=2000m/s┃┣━━━━━━━━┻━━━━━━━━┻━━━━━━━━┫┃第三层H=1400mV=2500m/s┃┗━━━━━━━━━━━━━━━━━━━━━━━━━━┛图1含油地质模型┏━━━━━━━━┓┃第一层H=1000mV=1500m/s┃┣━━━━━━━━┫┃V=2000m/sV=2500m/s┃┗━━━━━━━━┛在该地质模型中，作者设计了一个含油层，并使该含油层包裹在砂岩环境中，各层的速度及深度参数如模型中所示。

此外，为了更加直观地显示地震勘探的效果，本文提供了一个参照图2对比地质模型模型，该模型中没有设计含油速第37卷（2009）第7期结果，从地下许多反射界面发生反射时形成的地震子波，振幅大小决定于反射界面反射系数的绝对值，极性的正负决定于反射系数的正负，到达时间的先后取决于界面深度和覆盖层的波速‘1]震子波波形用S（t）表示，反射系数是双程垂直反射旅行时f的函数，用R（t）表示，地震记录厂（z）形成的物理过程在数学上就可表示为：

／’（f）一S（f）*R（f）一S（r）R（t-r）dr2）为了模型建立和计算过程中突出理论数值，去除了一些干扰因素，对一切衰减、噪声都不进行4）地震子波以平面波形式垂直入射到界面，各薄层的反射子波与地震子波形状相同，只是振幅及极性不同。

基于以上假设条件进行地震记录合成就必须已知地震子波以及地层的反射系数，而反射系数又主要由地层的波阻抗反映，所以必须首先获取地层的速度和密度资料。

速度资料可通过连续速度测井获得，密度资料可从密度测井获得，得不到密度资料时，可近似假定密度不变，以速度曲线代替波阻抗曲线来计算反射系数。

加德纳根据实际资料提出了一个由速度推算密度的经验公式：

10=0.23V'25（速度单位：

ft/s）告R~:

其中，V已知地震子波是合成地震记录的一个很重要的前提条件。

在已有的地震记录上选取地震子波的具体方法如下：

在地震记录上识别出单波，做出单波波形，再反复试验，检查找出符合实际的子波。

2）根据已总结出的地震子波特点，采取具有特殊数学表达式的波形表示，如雷克子波等。

3）使用非炸药震源时记录的震源子波波形。

4）利用实际地震记录，采取数字处理方法在一定的假设条件下求取地震子波。

6）在拥有声波测井资料和井旁地震记录X（t）时，反射系数曲线R（t），地震子波S（t），可由X（厂）=S（）*R（f）得子波的谱S（f）=X（厂）/R（／），再对上式作反傅氏变换得地震子波波形S（t）。

地质及地球物理学研究表明，地表之下的地质结构是极其复杂的。

地质构造是指地壳中的地层在岩层地壳运动的作用下发生变形与变位而遗留下来的形态…状均匀介质模型，即似没有发育在三个不同时期的水平层状沉积岩层。

同时，为了体现地震勘探在石油、天然气勘探领域的应用效果，作者结合实际地质资料模拟了地下储层含有石油的情况。

一般情况下，随着埋藏深度越来越大，岩石的地震波传播速度会越来越大。

深层的岩石地震波传播速度通常要大于浅层岩石的传播速度。

而流体的传播速度较岩石的传播速度要小很多。

比如，砂岩的地震波传播速度一般在1800m/s～4000m/s，石油的传播速度一般为1300m/s～1400m/s。

在模型建立的过程中，作者参照了实际地质构造及各种岩石的速度资料，确定了三层水平层状均匀介质的速度及深度参数。

储层含油地质模型如图1所示：

┏━┓V=1500m/s┣┳┫V=2000V=1300┻V=2500┗┛图1含油地质模型V=1500在该地质模型中，作者设计了一个含油层，并使该含油层包裹在砂岩环境中，各层的速度及深度参数如模型中所示。

此外，为了更加直观地显示地震勘探的效果，本文提供了一个参照134贾跃玮等：

基于MATLA13的地震正演模型实现第37卷度突变层，在地震模拟记录参数不边的条件下用来比照含油模型的效果。

模型构造及参数如图2。

4地震模型正演根据褶积理论，结合地质模型，作者在MAT-I。

AB环境中编写程序实现了二维地震正演。

首先，为了最大可能的与实际情况相符合，地震模型中使用的子波是稳定可实现的子波‘8]，如图3。

图3子波该子波是最小相位子波，有时称为前载子波，其能量集中在整个波形的前端。

由于大多数脉冲地震震源（如炸药震源）产生的原始脉冲是接近最小相位的，因此在地震正演模型中的地震子波选取一般都选择最小相位型子波‘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该模型为各层平行，无倾角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（2'pi*f*tt）.*exp（-b*tt）；视速度hvaO—h（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）;fori=2:

nlhva0（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）;endhva0—sqrt（hva0./hvbo）;hval=sqrt（hval./hvbl）;hva2一sqrt（hva2./hvb2）;hva0

（1）一v0

（1）：

hval（l）=vl（l）：

hva2

（1）一v2

（1）；以上程序段使用循环控制各个层速度变化fori=l:

nl设定反射面深度forj=2:

mh（i,j）=h（i,j-l）endendmk=round（m/2）；炮检距hvbO（mk）=0;hvbl（mk）=0;hvb2（mk）=o：

hvbl=（l:

m）*dx;hvbl=（l:

m）*dx;hvb2=（l:

m）*dx;fori=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；划分网隔fori=1：

nl一次反射波forj=1:

80ifj<=30;z=2*h（i,j）/hvaO（i）；t=round（z/dt）；x（t，j）一x（t，j）+rO（i）；elseifj>30&＆。

j<一50；z=2*h（i,j）/hval（i）；t—round（z/dt）：

x（t，j）一x（t，j）+r1（i）；elsez=2*h（i,j）/hva2（i）;t=round（z/dt）：

x（t，j）一x（t，j）+r2（i）；end卷度突变层，在地震模拟记录参数不边的条件下用来比照含油模型的效果。

模型构造及参数如图2。

I环境中编写程序实现了二维地震正演。

首先，3子波该子波是最小相位子波，有时称为前载子波，其能量集中在整个波形的前端。

由于大多数脉冲地震震源（如炸药震源）产生的原始脉冲是接近最小相位的，因此在地震正演模型中的地震子波选取一般都选择最小相位型子波‘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-60．b=30;子波参数wb=a0*sin（2'pi*f*tt）.*exp（-b*tt）视速度hvaO—h（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）;endhva0—sqrt（hva0./hvbo）;hval=sqrt（hval./hvbl）;hva2一sqrt（hva2./hvb2）;hva0

（1）一v0

（1）：

hval（l）=vl（l）forj=2:

mh（i,j）=h（