CN110161561B - 一种油气储层中的可控层位分阶层间多次波模拟方法 - Google Patents
一种油气储层中的可控层位分阶层间多次波模拟方法 Download PDFInfo
- Publication number
- CN110161561B CN110161561B CN201910468937.4A CN201910468937A CN110161561B CN 110161561 B CN110161561 B CN 110161561B CN 201910468937 A CN201910468937 A CN 201910468937A CN 110161561 B CN110161561 B CN 110161561B
- Authority
- CN
- China
- Prior art keywords
- horizon
- matrix
- wave field
- continuation
- wavefield
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 118
- 238000004088 simulation Methods 0.000 title claims abstract description 92
- 230000002441 reversible effect Effects 0.000 claims abstract description 30
- 230000003044 adaptive effect Effects 0.000 claims abstract description 7
- 125000004122 cyclic group Chemical group 0.000 claims abstract description 6
- 239000011159 matrix material Substances 0.000 claims description 65
- 239000010410 layer Substances 0.000 claims description 64
- 239000011229 interlayer Substances 0.000 claims description 24
- 239000013598 vector Substances 0.000 claims description 18
- 230000015572 biosynthetic process Effects 0.000 claims description 17
- 238000004364 calculation method Methods 0.000 claims description 14
- 238000005070 sampling Methods 0.000 claims description 12
- 238000001514 detection method Methods 0.000 claims description 5
- 230000005540 biological transmission Effects 0.000 claims description 4
- 238000003825 pressing Methods 0.000 claims description 4
- 230000008859 change Effects 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 238000012935 Averaging Methods 0.000 claims 2
- 238000010276 construction Methods 0.000 claims 2
- 230000008569 process Effects 0.000 abstract description 20
- 238000005516 engineering process Methods 0.000 abstract description 6
- 230000001629 suppression Effects 0.000 abstract description 5
- 238000002360 preparation method Methods 0.000 abstract 1
- 238000010586 diagram Methods 0.000 description 17
- 238000012937 correction Methods 0.000 description 12
- 230000000694 effects Effects 0.000 description 9
- 230000001174 ascending effect Effects 0.000 description 5
- 238000012545 processing Methods 0.000 description 5
- BVKZGUZCCUSVTD-UHFFFAOYSA-L Carbonate Chemical compound [O-]C([O-])=O BVKZGUZCCUSVTD-UHFFFAOYSA-L 0.000 description 3
- 235000019738 Limestone Nutrition 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 3
- 230000007547 defect Effects 0.000 description 3
- 239000006028 limestone Substances 0.000 description 3
- 230000000737 periodic effect Effects 0.000 description 3
- 206010013883 Dwarfism Diseases 0.000 description 2
- 230000002238 attenuated effect Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 125000004432 carbon atom Chemical group C* 0.000 description 2
- 239000003245 coal Substances 0.000 description 2
- 230000002349 favourable effect Effects 0.000 description 2
- 239000011435 rock Substances 0.000 description 2
- 238000000926 separation method Methods 0.000 description 2
- 238000004613 tight binding model Methods 0.000 description 2
- 239000004215 Carbon black (E152) Substances 0.000 description 1
- 241000282345 Meles Species 0.000 description 1
- 208000037386 Typhoid Diseases 0.000 description 1
- 230000002146 bilateral effect Effects 0.000 description 1
- 230000006866 deterioration Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 230000032669 eclosion Effects 0.000 description 1
- 229930195733 hydrocarbon Natural products 0.000 description 1
- 150000002430 hydrocarbons Chemical class 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 230000000670 limiting effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000002829 reductive effect Effects 0.000 description 1
- 238000001028 reflection method Methods 0.000 description 1
- 238000002310 reflectometry Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 201000008297 typhoid fever Diseases 0.000 description 1
- 238000010200 validation analysis Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公布了一种油气储层中的可控层位分阶层间多次波模拟方法，属于地震波场正演模拟领域，包括自适应延拓步长计算、二维反周期延拓边界反射处理和波场模拟中的双重层位约束。该方法基于自适应变步长波场延拓，以递归循环的方式实现分阶层间多次波的模拟；包括输入数据准备过程、递归循环波场模拟过程和压制边界反射过程。通过对模型添加双重层位约束，可以模拟指定地层产生的各阶层间多次波；利用二维反周期延拓方法压制波场延拓过程中产生的人工边界反射，得到了比传统衰减边界更佳的边界反射压制效果；提出自适应变步长波场延拓技术，大大提升了波场模拟的效率。
Description
技术领域
本发明属于地震波场正演模拟领域，涉及自适应变步长波场延拓、二维反周期延拓边界反射处理和波场模拟中的双重层位约束，尤其涉及一种油气储层中的可控层位分阶层间多次波模拟方法，是一种在炮集域模拟指定地层之间的一次波和各阶层间多次波波场的数值模拟方法。
背景技术
中国油气储层深度内普遍存在低速层，例如四川盆地的二叠系和奥陶系的低速泥岩页岩层，塔里木盆地的侏罗系煤层，伊利盆地的侏罗系煤层等。这些低速层的存在导致地震资料中层间多次波发育，严重影响对深部储层的处理和解释。有效识别和处理深部储层上覆地层产生的层间多次波是提高深部储层解释准确度必不可少的环节。
目前对层间多次波的识别方法主要可以分为两类，一类是数据驱动的层间多次波预测方法，另一类是模型驱动的波场模拟方法。一种数据驱动的预测类方法是由Ikelle在2006年提出[1]，后来由刘嘉辉等在2018年改进的虚同相轴方法[2]，该方法由拾取的一次波同相轴出发，预测与产生该一次波同相轴的层位相关的层间多次波，虚同相轴方法可以预测比较准确的层间多次波到时，但是不能预测准确的层间多次波振幅，且只能预测与某一层位相关的总的层间多次波。另一种数据驱动的预测类方法是基于Marchenko自聚焦的层间多次波预测，Meles等2014年首先将Marchenko方法引入到声波层间多次波的预测中[3]，Filho等2016年又将该方法推广至弹性波情况[4]，并在2018年将Marchenko方法应用于北海的实际数据中[5]，此方法首先利用Marchenko自聚焦恢复地下虚源点的上行格林函数波场和下行格林函数波场，然后结合干涉法将其用于构建层间多次波，这种方法计算量大，且也无法分阶预测层间多次波。
有限差分数值模拟是目前最常用的模型驱动类波场模拟方法，该方法具有较高的计算效率和精度[6]，工业上将其作为主要的波场模拟手段，但是该方法只能模拟总波场，无法将一次波场和层间多次波场分离。而层间多次波模拟方法可以实现一次波和层间多次波的分离，且模拟波场的时间和振幅都具有较高的精度。Kennett等在1974提出了一种基于反射率法的层间多次波模拟方法[7]，但该方法目前无法很好地适应介质的横向变化，且无法实现指定地层之间的层间多次波分阶模拟；Covey等1989年提出了一种基于射线追踪的层间多次波模拟方法[8]，该方法运算效率较高，能处理简单介质中的层间多次波模拟问题，但无法应对复杂变化的模型。最近，Berkhout在2014提出了一种基于波场延拓的全波场模拟方法[9]，弥补了有限差分方法只能一次性模拟所有波场的缺陷，可以分阶次模拟波场，且能够处理相对复杂的地质模型，Davydenko等2018年将其应用于实际数据全波场成像取得了一定的效果[10]，但是，该方法目前只能一次性模拟所有地层的各阶层间多次波，不能直接模拟给定地层之间各阶层间多次波，且计算效率有待提升。
参考文献：
[1]Ikelle,L.T..A construct of internal mutiples form surface dataonly:the concept of virtual seismic events.Geophysical JournalInternational.2006,164(2),383–393.
[2]Liu J.H.,Hu T.Y.,Peng G.X.,et al.Removal of internal multiples byiterative construction of virtual primaries.Geophysical JournalInternational.2018,215(1):81-101.
[3]Meles G.A.,Loer K.,Ravasi M.,et al.Internal multiple predictionand removal using Marchenko autofocusing and seismicinterferometry.Geophysics.2014,80(1):A7-A11.
[4]Filho C.A.D.C.,Meles G.A.,and Curtis A..Elastic internal multipleprediction using Marchenko and interferometric methods.SEG Technical ProgramExpanded Abstracts.2016,4545-4549.
[5]Filho C.A.D.C.,and Curtis A..Marchenko and interferometry basedmultiple attenuation of a North Sea field dataset.SEG Technical ProgramExpanded Abstracts.2018,4518-4522.
[6]Zhou H.,Liu Y.,and Wang J..Finite difference modeling withadaptive variable-length temporal and spatial operators.SEG Technical ProgramExpanded Abstracts.2018,4015-4019.
[7]Kennett B.L.N..Theoretical reflection seismograms for elasticmedia.Geophysical Prospecting.1979,27(2):301-321.
[8]Covey J.D.,Hron F.,and Peacock K.L..On the role of partial rayexpansion in the computation of ray synthetic seismograms for layeredstructures.Geophysical Prospecting.1989,37(8):907-923.
[9]Berkhout A.J..Review paper:An outlook on the future of seismicimaging,Part I:forward and reverse modeling.Geophysical Prospecting.2014,64(5):911-930.
[10]Davydenko M.,and Verschuur D.J..Including and using internalmultiples in closed-loop imaging-Field data examples.Geophysics.2018,83(4):R297-R305.
发明内容
针对上述现有技术存在的问题，本发明提供一种油气储层中的可控层位分阶层间多次波模拟方法，输入数据包括地质模型(本说明书中表示速度模型和密度模型的统称)和在地质模型中拾取的约束层位数据，通过构建层位约束矩阵、自适应延拓步长计算、震源设置、递归循环波场模拟和边界反射压制，实现指定地层之间的各阶层间多次波模拟，能够提升层间多次波波场模拟的效率，有利于辅助识别层间多次波。
本发明方法基于全波场模拟方法，在算法中添加双重层位约束，即指定产生下行反射的地层和产生上行反射的地层，使层间多次波的模拟更加灵活；引入了时间-空间域二维反周期延拓方法处理波场模拟过程中的人工边界反射；提出了自适应变步长波场延拓技术，用以提升层间多次波波场模拟的效率。并将本方法应用到中国某油田地区的实际模型中，对实际数据中的层间多次波起到了很好的辅助识别作用。
本发明的核心是：基于波场延拓全波场模拟技术，在算法中添加双重层位约束指定上行反射和下行反射的地层，更加自由灵活地模拟指定主要地层之间的层间多次波；引入二维反周期延拓方法压制波场模拟过程中产生的人工边界反射，压制人工边界反射对模拟有效信号的干扰；提出了自适应变步长波场延拓技术，根据地质模型的在纵向上的复杂度，在不复杂的地层内使用大步长，在复杂的地层内使用小步长，在不影响模拟精度的情况下提升计算效率。相比于常规的有限差分方法，本发明能在保证模拟精度的情况下模拟指定地层之间的各阶层间多次波波场，而不是仅仅模拟一个总的波场，有利于对数据中层间多次波进行更加精细的识别和处理。
本发明提供的技术方案如下：
一种油气储层中的可控层位分阶层间多次波模拟方法，输入数据包括地质模型和在地质模型中拾取的约束层位数据，通过构建层位约束矩阵、自适应延拓步长计算、震源设置、递归循环波场模拟和边界反射压制，实现指定地层之间的各阶层间多次波模拟；包括以下步骤(图1)：
A.数据准备：
A1.构建层位约束矩阵：根据实际需要，在地质模型中拾取发生下行反射的层位和发生上行反射的层位，用于构建下行反射层位约束矩阵
和上行反射层位约束矩阵
A2.自适应延拓步长计算：根据地质模型的复杂度，计算地质模型随深度变化的延拓步长；
A3.震源设置：根据时间采样率、时间记录长度、检波点数、震源位置及震源子波定义原始震源矩阵s，对s进行二维反周期延拓形成反周期延拓后的反周期延拓震源矩阵
B.递归循环波场模拟：
B1.初始化，将原始震源向量
赋值给地表的下行波场向量
并将地下各个层位下方的上行波场向量
赋值为零向量；
B2.将地表的下行波场向量
根据A2中计算的自适应延拓步长逐步向下延拓至模型底部，得到各个层位上方的下行波场向量
和层位下方的下行波场向量
及模型底部层位上方的下行波场
B3.将B2中得到的模型底部上方的下行波场
乘以模型底部层位的下行反射系数矩阵
得到上行波场
B4.将模型底部层位上方的上行波场
逐步向上延拓至地表，得到各个层位上方的上行波场向量
和层位下方的上行波场向量
及地表层位下方的上行波场为
即我们模拟的采集数据；
B5.递归执行B2-B4，每循环一次得到更高一阶的层间多次波波场，直到得到所需阶次的层间多次波波场；
B6.用反周期延拓震源
替换原始震源
执行B2-B5，模拟得到二维反周期延拓后的层间多次波波场；
C.压制边界反射：由于傅里叶变换的周期性假设，将B5和B6中对应阶次的层间多次波波场求平均，得到压制边界反射后的层间多次波波场。
本发明的有益效果是：
传统有限差分数值模拟方法通过数值求解波动方程模拟地震波经过介质传播后的总波场，包括所有一次波和多次波，无法得到阶次分离的波场。而本发明基于自适应变步长波场延拓，以递归循环的方式实现一次波场和层间多次波波场的模拟。本发明的优点在于：
(一)通过波场延拓，以递归循环的方式实现层间多次波的分阶模拟，克服了常规有限差分方法只能模拟总波场的缺点，可以分阶次模拟一次波波场和层间多次波波场；
(二)添加双重层位约束，实现模拟指定层位之间的层间多次波，有利于层间多次波的精细识别；
(三)引入二维反周期延拓，达到了比常规衰减边界更佳的人工边界反射压制效果；
(四)提出自适应变步长延拓技术，减少完成波场模拟所需要的延拓步数，提升波场模拟效率。
附图说明
图1是本发明提供的一种油气储层中可控层位分阶层间多次波模拟方法的流程框图。
图2是波场延拓步长划分示意图；
其中，虚线表示每一步波场延拓的层位；(a)为固定步长波场延拓步数划分示意图；(b)为自适应变步长波场延拓步数划分示意图。
图3是波场在传播过程中的关系示意图；
其中，(a)为下行延拓过程中的波场关系示意图；(b)为上行延拓过程中的波场关系示意图。
图4是层间多次波模拟过程中的波场循环示意图；
其中，黑色倒三角表示震源，黑色圆点表示激发的二次源；(a)为第一次下行循环示意图；(b)为第一次上行循环示意图；(c)为第二次下行循环示意图；(d)为第二次上行循环示意图。
图5是本发明方法验证中用到的四层速度模型。
图6是本发明实施例中对比人工边界压制效果用到的炮集数据图；
其中，(a)为原始未压制边界反射单炮道集；(b)为常规衰减边界压制边界反射后道集；(c)为二维反周期延拓压制边界反射后道集；(a)中的黑色箭头指示人工边界反射产生的干扰，(b)中的黑色箭头指示传统衰减边界方法压制边界反射的残余。
图7是本发明实施例中对比本发明方法模拟结果和常规有限差分方法结果的炮集数据图；
其中，(a)为有限差分方法模拟的波场；(b)为本发明模拟的总波场；(c)为本发明模拟的一阶层间多次波波场。
图8是本发明实施例中对比本发明方法模拟结果和常规有限差分方法结果的单道数据图；
其中，(a)为比较本发明模拟的一次波和有限差分结果图；(b)为比较本发明模拟的一阶层间多次波和有限差分结果图；(c)为比较本方法模拟的二阶层间多次波和有限差分结果图。
图9是本发明在实施例中用到的某油田地区的二维速度模型图；
其中，(a)为完整速度模型图；(b)为完整速度模型图中虚线方框范围的放大；(b)中虚线标注三组产生层间多次波的主要地层组，分别用A、B、C标注，数字1、2、3所标示的传播路径表示A、B、C三组地层之间产生层间多次波的三种不同路径，路径中的粗线表示波场在地层组中的反射是该地层组的总反射效应。
图10是本发明在实施例中用某油田地区的二维速度模型模拟的数据和实际采集数据第940个CMP道集数据的对比图；
其中，(a)为动校正前实际观测CMP道集；(b)为动校正前模拟总波场CMP道集；(c)为动校正前模拟一阶层间多次波CMP道集；(d)为动校正后实际观测CMP道集；(e)为动校正后模拟总波场CMP道集；(f)为动校正后模拟一阶层间多次波CMP道集；(c)图中的黑色箭头标注主要层间多次波的位置，(d)图和(f)图中白色箭头和灰色箭头标注短程层间多次波，黑色箭头标注长程层间多次波。
图11是本发明在实施例中用某油田地区的二维速度模型模拟添加层位约束后的第940个CMP道集数据图；
其中，(a)为未添加层位约束前模拟的总波场CMP道集；(b)为不允许层位组B发生上行反射和下行反射后模拟总波场CMP道集；(c)为(a)和(b)的CMP差异道集；(d)为仅允许层位组B发生上行反射和下行反射模拟得到的一阶层间多次波CMP道集；(e)为允许层位组B发生上行反射层位组A发生下行反射模拟得到的一阶层间多次波CMP道集；(f)为允许层位B发生上行反射和下行反射，并且允许层位组C发生上行反射模拟得到的一阶层间多次波CMP道集；图(a)和(b)中的箭头指示层间多次波的主要能量；图(c)中白色箭头指示允许层位组B发生上行反射和下行反射模拟得到的一阶层间多次波，黑色箭头指示允许层位组B发生上行反射层位组A发生下行反射模拟得到的一阶层间多次波，灰色箭头指示允许层位组B发生上行反射和下行反射，且允许层位组C发生上行反射模拟得到的一阶层间多次波。
图12是本发明在实施例中用某油田地区的二维速度模型模拟的一阶层间多次波叠加剖面和实际数据叠加剖面的对比；
其中，(a)为实际地震数据叠加剖面图；(b)为本发明方法模拟的一阶层间多次波叠加剖面图；图中白色箭头指示仅允许层位组B发生上行反射和下行反射模拟得到的一阶层间多次波同相轴，黑色箭头指示允许层位组B发生上行反射层位组A发生下行反射模拟得到的一阶层间多次波同相轴，灰色箭头指示允许层位组B发生上行反射和下行反射，且允许层位组C发生下行反射模拟得到的一阶层间多次波同相轴。
具体实施方式
下面结合附图，通过实施例进一步描述本发明，但不以任何方式限制本发明的范围。
本发明基于全波场模拟方法，在算法中添加双重层位约束，即指定产生下行反射的地层和产生上行反射的地层，使层间多次波的模拟更加灵活；引入了时间-空间域二维反周期延拓方法处理波场模拟过程中的人工边界反射；提出了自适应变步长波场延拓技术，用以提升层间多次波波场模拟的效率。并将本发明方法应用到中国某油田地区的实际模型中，对实际数据中的层间多次波起到了很好的辅助识别作用。本发明输入数据包括地质模型和在地质模型中拾取的约束层位数据，通过构建层位约束矩阵、自适应延拓步长计算、震源设置、递归循环波场模拟和边界反射压制，实现指定地层之间的各阶层间多次波模拟。
图1是本发明提供的一种可控层位分阶层间多次波模拟方法流程框图，方法包含以下步骤：
A.数据准备：
A1.构建层位约束矩阵，包括下行反射层位约束矩阵
和上行反射层位约束矩阵
根据实际地震数据或者测井数据，估计可能的产生层间多次波主要能量的上行反射层位和下行反射层位，在地质模型中拾取发生下行反射的层位和发生上行反射的层位。设地质模型的纵向网格点数为Nz,横向网格点数为Nx,则构建两个大小为Nz×Nx的矩阵，分别作为上行约束模型矩阵Gu和下行约束模型矩阵Gd。Gd和Gu的计算方式为：若Gd的第n行第m列位于地质模型中拾取的下行反射地层上，则(Gd)n,m＝1，否则(Gd)n,m＝0；若Gu的第n行第m列位于地质模型中拾取的下行反射地层上，则(Gu)n,m＝1，否则(Gu)n,m＝0。
上式中，
为下行反射层位约束矩阵；
为上行反射层位约束矩阵，上标n为层位位置。本说明书中矩阵括号外的第一个下标表示矩阵的行号，第二个下标表示矩阵的列号，此处行号i的取值范围为1到Nx的整数，行号n的取值范围为1到Nz的整数，列号i的取值范围也是1到Nx的整数。
A2.计算得到自适应延拓步长：根据地质模型的复杂度，对于构造简单的地层采用大步长，对于构造复杂的地层采用小步长，计算得到模型随地层深度变化的延拓步长。
对于同一个地质模型，基于波场延拓的层间多次波模拟的计算效率在一定程度上取决于延拓步长的大小。图2是一个2000m*1600m的地质模型，水平方向和深度方向的网格间距都是50m，常规的基于波场延拓的模拟方法采用固定步长，即延拓步长为50m不变，其延拓步长划分如图2(a)中的虚线所示，本例中一次循环需要的延拓步数为33步；在不降低模拟精度的前提下，可以采取自适应变步长波场延拓的方式提高波场模拟的效率，根据地下介质的复杂程度，在构造简单的地方采用大步长，在构造复杂的地方采用小步长。
上式中，Cn为地质模型第n个层位的变步长因子，
和
分别为地质模型第n个层位和第n-1个层位的速度向量。当Cn＝0时，将第n个层位与第n-1个层位合并为一个延拓步长，当Cn≠0时，将第n个层位单独作为一个延拓步长。对图2所示的模型做自适应步长划分后的结果如图2(b)所示，本例中一次循环需要的延拓步数为14步。每一步波场延拓的计算效率与步长无关，则在整个波场模拟的过程中，所需的步数与计算量成线性负相关的关系，即所需步数越少，总的计算量越小，计算效率越高。图2为了显示方便，设计的模型网格间距为50m，在实际情况下，当网格间距越小时，自适应变步长方法对模拟效率的提升越显著。
A3.震源设置：根据时间采样率、时间记录长度、检波点数、震源位置及震源子波定义原始震源矩阵s，对s进行二维反周期延拓，形成反周期延拓后的反周期延拓震源矩阵
具体计算过程如下：
设时间采样率为dt，时间记录长度为T，计算得到时间采样点数
设检波点数为Nr，则原始震源矩阵为一个大小Nt×Nr的矩阵s；设震源位置为x，地震子波为雷克子波，子波主频为fp，考虑到雷克子波的零相位特性，设子波长度为Nw为奇数,则原始震源矩阵s的第x列第i行的值为：
上式中等式左边的下标i和x分别表示原始震源矩阵s的行号和列号。根据原始震源矩阵s，可以得到K阶二维反周期延拓后的反周期延拓震源矩阵
具体计算公式为：
上式中k1和k2满足0≤k1,k2<K，k1,k2,K均为非负整数，通常情况下取K＝2便可得到比较理想的边界压制效果，I表示水平方向空间采样点数，J表示时间采样点数，上式中所有自变量均为整数。将原始震源矩阵s变换到频率域，然后取其单频向量作为步骤B1中的输入震源
将反周期延拓震源矩阵
变换到频率域，然后取其单频向量作为步骤B6中的替换输入震源
B.递归循环波场模拟：
本发明中，以层位n表示地质模型纵向第n层网格。将层间多次波模拟过程中层位n处的波场分为四类，即层位上方的下行波场和上行波场，层位下方的下行波场和上行波场，在频率域分别用
单频列向量表示，向量中的元素对应不同的水平坐标，P表示层位上方波场，Q表示层位下方波场，上标n表示层位位置，下标u表示上行波，下标d表示下行波。波场在传播过程中的关系如图3所示，图3(a)表示下行延拓过程中的波场关系，图3(b)表示上行延拓过程中的波场关系。
在下行延拓的过程中，层位n-1下方的下行波场
经过下行延拓至层位n上方，得到
层位n上方的下行波场
经过层位n的透射波场加上层位n下方的上行波场
经层位n的反射波场得到层位n下方的下行波场
在上行延拓的过程中，层位n+1上方的上行波场
经过上行延拓至层位n下方，得到
层位n下方的上行波场
经过层位n的透射波场加上层位n上方的下行波场
经层位n的反射波场得到层位n上方的上行波场
上式中，
表示下行延拓矩阵算子，将层位n-1处的波场延拓至层位n，
表示上行延拓矩阵算子，将层位n+1处的波场延拓至层位n。
表示层位n处的下行透射系数矩阵，
表示层位n处的上行反射系数矩阵，
表示层位n处的上行透射系数矩阵，
表示层位n处的下行反射系数矩阵，
为下行反射层位约束矩阵，
上行反射层位约束矩阵。
B.递归循环波场模拟，包括如下步骤：
B1：初始化，
n>0时
B2：将
代入(式6)中得到
再将
和
代入(式7)中得到
将得到的波场循环代入(式6)和(式7)，直到将波场延拓至地质模型底部，得到所有层位上方和下方的下行波场
和
及模型底部层位上方的下行波场
此过程如图4(a)所示；
本发明中用层位e表示模型最深部的层位。
B3：将B2中得到的模型底部上方的下行波场
乘以模型底部层位的反射系数矩阵
得到上行波场
B4：将模型底部层位上方的上行波场
代入(式8)得到层位e-1下方的上行波场
将B2中得到的
和前面得到的
代入(式9)得到层位e-1上方的上行波场
将得到的波场循环代入(式8)和(式9)，直到将波场延拓至地表，得到所有层位上方和下方的上行波场
和
及地表的上行波场
即采集数据，此过程如图4(b)所示；
B5：递归执行B2-B4，每循环一次得到更高一阶的层间多次波波场，直到得到所需阶次的层间多次波波场，此过程如图4(c)和4(d)所示。
B6.用反周期延拓震源
替换原始震源
执行B1-B5，模拟得到二维反周期延拓后的层间多次波波场；
C.压制边界反射：由于傅里叶变换的周期性假设，将B5和B6中对应阶次的层间多次波波场求平均，得到压制边界反射后的层间多次波波场。
本发明在一个简单的四层模型上对比了本发明方法所采用的二维反周期延拓方法和常规的衰减边界压制人工边界反射方法的效果，并比较采用本方法方法模拟得到的波场和采用有限差分方法得到的波场的一致性。地质模型水平长度为3000m，深度为1750m，其速度模型如图5所示。
在模拟过程中分别采用常规的衰减边界和本发明提出的二维反周期延拓方法压制人工边界反射，模拟结果如图6所示，图6(a)为原始不压制人工边界反射的单炮道集，图6(b)为常规衰减边界压制人工边界反射后的单炮道集，图6(c)为本发明采用的二维反周期延拓法压制人工边界反射后的单炮道集。
由于波场延拓过程中正反傅里叶变换的周期性假设，导致波场在模型边界发生边界反射并传导到时间轴，使得边界反射严重干扰反射波场，如图6(a)中箭头所示，使得模拟结果让人无法接受；常规的衰减边界是在波场模拟的过程中对衰减边界内的波场乘以一个e指数衰减项，以达到衰减边界反射的目的，本发明添加衰减边界压制边界反射的效果如图6(b)所示，与图6(a)相比，衰减边界压制了大部分的边界反射，但是仍然有部分的边界反射残留，如图6(b)中箭头所示；而对比图6(b)和6(c)可以看到，本发明的二维反周期延拓方法比常规的衰减边界方法有更佳的压制边界反射的效果。
为了进一步验证本方法的正确性，将本发明方法模拟得到的波场和有限差分模拟的结果进行比较如图7所示。图7(a)为有限差分模拟的波场，图7(b)为本发明方法模拟的总波场，图7(c)为本发明方法模拟的一阶层间多次波场。
对比图7(a)和7(b)可以看到，本发明方法模拟得到的波场与有限差分方法模拟的波场具有很好的一致性，同时本发明方法还克服了有限差分方法无法分阶模拟波场的缺陷。为了更清晰地比较本发明方法和有限差分方法的结果，我们在零偏移距单道上对比本发明方法的结果和有限差分方法的结果，比较结果如图8所示，图8(a)比较本发明方法模拟的一次波和有限差分结果，图8(b)比较本发明方法模拟的一阶层间多次波和有限差分结果，图8(c)比较本发明方法模拟的二阶层间多次波和有限差分结果，其中实线表示有限差分模拟结果，虚线为本发明方法的模拟结果。
图8中P表示一次波，P后面的数字表示产生一次波的层位，M表示层间多次波，M后面的数字表示发生下行反射的层位，一个数字表示仅发生一次下行反射，即一阶层间多次波，两个数字表示分别在两个数字所指的层位发生下行反射，即二阶层间多次波。从图8的对比结果可以看到，采用本发明的方法模拟得到的一次波和各阶层间多次波均与有限差分模拟结果很好的一致性，从而验证了本方法的正确性。
为了证明本方法在实际地质模型中的有效性，我们将本发明的方法应用到某油田地区的二维模型上，二维地质模型通过对该地区的实际地震数据和测井数据建模得到。该区目标层位于深度约5.5km-6.5km的奥陶系碳酸盐岩储层，目标区域上覆石炭系灰岩地层，目标层下伏寒武系基底，石炭系灰岩地层之上发育广泛的侏罗系低速层，侏罗系低速层与石炭系灰岩层之间存在显著的波阻抗差，产生能量较强的层间多次波，对奥陶系碳酸盐岩储层的反射波信号形成严重的干扰，正确识别处理干扰奥陶系碳酸盐储层反射信号的层间多次波能量有助于提升目标层的解释准确度。该区二维速度模型如图9所示，图9(a)为整体速度模型，图9(b)为9(a)中虚线方框的放大模型，图9(b)中字母A、B、C标注三组地层，如图9(b)中虚线所示，数字1、2、3所标示的传播路径表示A、B、C三组地层之间产生层间多次波的三种不同路径，路径中的粗线表示波场在地层组中的反射是该地层组的总反射效应。
模拟数据最大偏移距为6175m，道间距50m，双边采集，模拟炮数为900炮，起始炮水平坐标位置为6175m，炮间距为50m，时间记录长度为7s，时间采样率为0.004s。实际数据第940个CMP道集和模拟数据相应CMP道集如图10所示，图10(a)-(c)分别为动校正前实际观测波场的CMP道集、模拟总波场数据CMP道集、模拟一阶层间多次波CMP道集，图10(d)-(f)分别为动校正后实际观测波场的CMP道集、模拟总波场CMP道集、模拟一阶层间多次波CMP道集，从图10(c)可以看到，层间多次波能量主要集中在2.8s-3.0s和3.2s-3.5s，如图10(c)中的黑色箭头所指；2.8s-3.0s的层间多次波动和3.4s-3.5s的层间多次波校时差与一次波相近，难以从原始CMP道集识别出来，如图10(d)和10(f)中的白色箭头和灰色箭头所指；3.2s-3.4s的层间多次波动校时差明显大于一次波动校时差，相比之下更易于从原始CMP道集识别出来，如图10(d)和10(f)中的黑色箭头所指。
为了分析影响储层反射信号的主要层间多次波能量的来源，我们添加双重层位约束控制波场模拟过程产生中产生上行反射和下行反射的层位，其结果如图11所示。图11(a)为未添加层位约束时模拟的总波场动校正后的CMP道集，可以看到数据中存在明显的层间多次波能量，如图11(a)中的箭头所示；图11(b)为添加约束不允许图9(b)中层位组B发生上行反射和下行反射模拟得到的总波场动校正后的CMP道集，对比图11(a)和11(b)可以看到，添加约束后，层间多次波的主要能量明显衰减了，由此可以判断，数据中的层间多次波主要能量都经过了层位组B的反射，图11(c)是图11(a)和11(b)CMP道集数据的差异。更进一步，我们仅允许层位组B发生上行反射和下行反射得到了图11(d)所示的动校正后的一阶层间多次波CMP道集，仅允许层位组A发生下行反射和层位组B发生上行反射，得到了图11(e)所示的动校正后的一阶层间多次波CMP道集，仅允许层位组B发生下行反射和上行反射，且允许层位组C发生上行反射，得到了图11(f)所示的动校正后的一阶层间多次波CMP道集，图11(d)(e)(f)中的层间多次波的射线路径分别如图9(b)中数字1、2、3所标注的路径。图11(d)所示的层间多次波由地层组B内的多次反射产生，图11(f)所示的层间多次波由地层组B内的多次反射波经地层组C上行反射后形成，这两组层间多次波可以归类为短程多次波，与一次波动校时差相近，在数据中往往难以识别。图11(e)所示的层间多次波由地层组A和地层组B之间的多次反射产生，这组层间多次波反射路径长，一次波与多次波动校时差差异较显著，可以归类为长程多次波，相比之下在数据中更容易识别。
为了分析层间多次波在叠加剖面上的特征，我们将本方法模拟得到的层间多次波叠加剖面和实际数据叠加剖面进行对比如图12所示，图12(a)为实际数据叠加剖面，图12(b)为本方法模拟的层间多次波叠加剖面。图中白色箭头所指为地层组B内的多次反射产生的层间多次波同相轴，黑色箭头所指为地层组A和地层组B之间的多次反射产生的层间多次波同相轴，灰色箭头所指为地层组B内的多次反射波经地层组C上行反射产生的层间多次波同相轴。层间多次波的存在导致叠加剖面的同相轴分辨率降低和连续性变差；图12(a)中白色箭头所指部位，一次波同相轴变宽并带有拖尾现象，分辨率明显下降；图12(a)中黑色箭头所指和灰色箭头所指同相轴难以追踪，连续性显著变差。可见，为了提高储层解释的准确度，避免假同相轴的干扰，必须进行有效的层间多次波识别和处理。
需要注意的是，公布实施例的目的在于帮助进一步理解本发明，但是本领域的技术人员可以理解：在不脱离本发明及所附权利要求的精神和范围内，各种替换和修改都是可能的。因此，本发明不应局限于实施例所公开的内容，本发明要求保护的范围以权利要求书界定的范围为准。
Claims (7)
1.一种油气储层中的可控层位分阶层间多次波模拟方法，通过构建层位约束矩阵、自适应延拓步长计算、震源设置、递归循环波场模拟和边界反射压制，实现油气储层中指定地层之间的各阶层间多次波模拟；包括以下步骤：
A.数据准备，执行如下操作：
A1.构建层位约束矩阵：
在地质模型中拾取发生下行反射的层位和发生上行反射的层位，用于构建下行反射层位约束矩阵
和上行反射层位约束矩阵
n为层位位置；d和u分别表示下行和上行；
A2.自适应延拓步长计算，即计算地质模型随深度变化的延拓步长；
A3.震源设置：
根据时间采样率、时间记录长度、检波点数、震源位置及震源子波定义原始震源矩阵s，对s进行二维反周期延拓，形成反周期延拓后的反周期延拓震源矩阵
B.递归循环波场模拟，执行如下操作：
B1.初始化：将原始震源
赋值给地表的下行波场
并将地下各个层位下方的上行波场
赋值为零向量；
B2.将地表的下行波场
根据步骤A2计算得到的自适应延拓步长逐步向下延拓至地质模型底部，得到各个层位上方的下行波场
和层位下方的下行波场
及模型底部层位上方的下行波场
层位e表示模型最深部的层位；
B3.将B2中得到的模型底部层位上方的下行波场
乘以模型底部层位的下行反射系数矩阵
得到上行波场
B4.将模型底部层位上方的上行波场
逐步向上延拓至地表，得到各个层位上方的上行波场
和层位下方的上行波场
及地表层位下方的上行波场为
即模拟的采集数据；
B5.递归执行B2-B4，每循环一次得到更高一阶的层间多次波波场，直到得到所需阶次的层间多次波波场；
B6.反周期延拓波场模拟，用反周期延拓震源
替换原始震源
执行B2-B5，模拟得到二维反周期延拓后的层间多次波波场；
C.压制边界反射：将B5和B6中对应阶次的层间多次波波场求平均，得到压制边界反射后的层间多次波波场；
通过步骤A～C，实现油气储层中的可控层位分阶层间多次波模拟。
2.如权利要求1所述的层间多次波模拟方法，其特征是，步骤A1所述构建层位约束矩阵，具体是根据实际地震数据或测井数据，估计可能的产生层间多次波主要能量的上行反射层位和下行反射层位，在地质模型中分别拾取估计的发生下行反射的层位和发生上行反射的层位；包括如下步骤：
设地质模型的纵向网格点数为Nz，横向网格点数为Nx，构建两个大小为Nz×Nx的矩阵，分别作为上行约束模型矩阵Gu和下行约束模型矩阵Gd；
Gd和Gu的构建方式为：若Gd的第n行第m列位于地质模型中拾取的下行反射地层上，则(Gd)n,m＝1，否则(Gd)n,m＝0；若Gu的第n行第m列位于地质模型中拾取的下行反射地层上，则(Gu)n,m＝1，否则(Gu)n,m＝0；
式中，
为下行反射层位约束矩阵；
为上行反射层位约束矩阵，上标n为层位位置；矩阵括号外的下标分别表示矩阵的行号和列号；m的取值范围为1到Nx的整数，n的取值范围为1到Nz的整数。
3.如权利要求2所述的层间多次波模拟方法，其特征是，步骤A2所述自适应延拓步长计算，具体通过式3计算得到：
式中，Cn为模型第n个层位的变步长因子，
和
分别为模型第n个层位和第n-1个层位的速度向量；当Cn＝0时，将第n个层位与第n-1个层位合并为一个延拓步长；当Cn≠0时，将第n个层位单独作为一个延拓步长。
4.如权利要求2所述的层间多次波模拟方法，其特征是，步骤A3的震源设置具体通过以下方式实现：
设时间采样率为dt，时间记录长度为T，计算得到时间采样点数
设检波点数为Nr，则原始震源矩阵为一个大小Nt×Nr的矩阵s；设震源位置为x，震源子波为雷克子波，子波主频为fp，设子波长度为Nw为奇数,则原始震源矩阵s的第x列第i行的值表示为式4：
式中，等式左边的下标i和x分别表示原始震源矩阵s的行号和列号；
根据原始震源矩阵s，通过式5计算得到K阶二维反周期延拓后的反周期延拓震源矩阵
式中，1≤i≤I,1≤j≤J；k1和k2满足0≤k1<K，0≤k2<K，K为常数；I表示水平方向空间采样点数，J表示时间采样点数；
将原始震源矩阵s变换到频率域，再取其单频向量作为步骤B1中的输入震源
将反周期延拓震源矩阵
变换到频率域，再取其单频向量作为步骤B6中的输入震源
5.如权利要求2所述的层间多次波模拟方法，其特征是，步骤B2的向下延拓递归循环中的每一步下行延拓具体通过如下方式实现：
在下行延拓的过程中，层位n-1下方的下行波场
经过下行延拓至层位n上方，得到
表示为式6：
层位n上方的下行波场
经过层位n的透射波场加上层位n下方的上行波场
经层位n的反射波场得到层位n下方的下行波场
其中，
表示下行延拓矩阵算子，将层位n-1处的波场延拓至层位n；
表示层位n处的下行透射系数矩阵，
表示层位n处的上行反射系数矩阵。
6.如权利要求2所述的层间多次波模拟方法，其特征是，步骤B4的向上延拓递归循环中的每一步上行延拓具体通过如下方式实现：
在上行延拓的过程中，层位n+1上方的上行波场
经过上行延拓至层位n下方，得到
表示为式8：
层位n下方的上行波场
经过层位n的透射波场加上层位n上方的下行波场
经层位n的反射波场得到层位n上方的上行波场
7.如权利要求1所述的层间多次波模拟方法，其特征是，步骤C的压制边界反射，具体通过截取原始波场和反周期延拓波场的重叠部分取平均实现。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910468937.4A CN110161561B (zh) | 2019-05-31 | 2019-05-31 | 一种油气储层中的可控层位分阶层间多次波模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910468937.4A CN110161561B (zh) | 2019-05-31 | 2019-05-31 | 一种油气储层中的可控层位分阶层间多次波模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110161561A CN110161561A (zh) | 2019-08-23 |
| CN110161561B true CN110161561B (zh) | 2020-09-08 |
Family
ID=67630830
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910468937.4A Active CN110161561B (zh) | 2019-05-31 | 2019-05-31 | 一种油气储层中的可控层位分阶层间多次波模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110161561B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111025386B (zh) * | 2019-12-13 | 2020-11-17 | 成都理工大学 | 一种无分离假象的纵横波分离方法 |
| CN115236730B (zh) * | 2022-06-22 | 2024-05-17 | 北京大学 | 一种层间多次波傅里叶有限差分的地震波场偏移成像方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102169189A (zh) * | 2011-01-19 | 2011-08-31 | 中国海洋石油总公司 | 深水层间多次波消除方法 |
| WO2012021218A2 (en) * | 2010-08-10 | 2012-02-16 | Geco Technology B.V. | Attenuating internal multiples from seismic data |
| CN109507722A (zh) * | 2017-09-15 | 2019-03-22 | 中国石油化工股份有限公司 | 基于模型及双波场延拓的层间多次波预测方法及系统 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8773949B2 (en) * | 2009-11-03 | 2014-07-08 | Westerngeco L.L.C. | Removing noise from a seismic measurement |
-
2019
- 2019-05-31 CN CN201910468937.4A patent/CN110161561B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2012021218A2 (en) * | 2010-08-10 | 2012-02-16 | Geco Technology B.V. | Attenuating internal multiples from seismic data |
| CN102169189A (zh) * | 2011-01-19 | 2011-08-31 | 中国海洋石油总公司 | 深水层间多次波消除方法 |
| CN109507722A (zh) * | 2017-09-15 | 2019-03-22 | 中国石油化工股份有限公司 | 基于模型及双波场延拓的层间多次波预测方法及系统 |
Non-Patent Citations (3)
| Title |
|---|
| Free-surface and internal multiple elimination in one step without adaptive subtraction;Lele Zhang et al.;《GEOPHYSICS》;20190228;第84卷(第1期);第A7-A11页 * |
| 基于反周期扩展边界方法的单程波正演模拟;沈铭成等;《中国西部科技》;20110531;第10卷(第14期);第45-46页 * |
| 层间多次波辨识与压制技术的突破及意义-以四川盆地GS1井区震旦系灯影组为例;甘利灯等;《石油勘探与开发》;20181231;第45卷(第6期);第960-971页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110161561A (zh) | 2019-08-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| AU2016331881B8 (en) | Q-compensated full wavefield inversion | |
| CN103238158B (zh) | 利用互相关目标函数进行的海洋拖缆数据同时源反演 | |
| CA2575274C (en) | System for attenuation of water bottom multiples in seismic data recorded by pressure sensors and particle motion sensors | |
| US5995905A (en) | Source signature determination and multiple reflection reduction | |
| EP3507626B1 (en) | Attenuation of multiple reflections | |
| US7257492B2 (en) | Handling of static corrections in multiple prediction | |
| CN108845351A (zh) | 一种vsp地震资料转换波全波形反演方法 | |
| US6999879B2 (en) | Method for controlling seismic coverage using decision theory | |
| CN111766628A (zh) | 一种预条件的时间域弹性介质多参数全波形反演方法 | |
| EP4031910B1 (en) | Noise attenuation methods applied during simultaneous source deblending and separation | |
| CN110161561B (zh) | 一种油气储层中的可控层位分阶层间多次波模拟方法 | |
| CN112946732A (zh) | 海上立体观测系统联合压制多次波单缆处理方法及系统 | |
| CN103119472A (zh) | 利用同时和顺序源方法进行全波形反演的混合方法 | |
| Sherwood et al. | Hybrid tomography based on beam migration | |
| CN112748463A (zh) | 一种基于深度学习照明分析的局部偏移成像方法 | |
| AU739128B2 (en) | A method of seismic processing, and in particular a 3D seismic prospection method implementing seismic data migration | |
| Kelamis et al. | Data-driven internal multiple attenuation-Applications and issues on land data | |
| Dutta et al. | Practical strategies for interbed multiple attenuation | |
| CN112147686B (zh) | 多期发育火成岩叠前深度偏移成像速度建模方法及系统 | |
| Cassell et al. | Interactive VSP-CDP mapping in complex media | |
| CN112946733A (zh) | 海上立体观测系统联合压制多次波多缆处理方法及系统 | |
| EP3938810B1 (en) | Signal recovery during simultaneous source deblending and separation | |
| CN113064205B (zh) | 菲涅尔带约束的浅水多次波衰减方法 | |
| Hu et al. | Addressing Complex Structural Imaging Challenges in Southern Caspian Sea Through Modern OBN Acquisition and Processing | |
| CN114428335A (zh) | 一种省时层间多次波预测方法和压制方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |