CN104133987B - 一种碳酸盐岩储层的逆时偏移方法 - Google Patents
一种碳酸盐岩储层的逆时偏移方法 Download PDFInfo
- Publication number
- CN104133987B CN104133987B CN201410325458.4A CN201410325458A CN104133987B CN 104133987 B CN104133987 B CN 104133987B CN 201410325458 A CN201410325458 A CN 201410325458A CN 104133987 B CN104133987 B CN 104133987B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msup
- mfrac
- msubsup
- pml
- 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.)
- Expired - Fee Related
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种碳酸盐岩储层的逆时偏移方法，该方法采用针对碳酸盐岩的缝洞孔隙构造常伴随有由于构造运动所形成的断裂和大角度的形态扭曲，运用逆时偏移方法对其进行精确刻画，本发明有限差分方法求解波动方程进行地震波场模拟，波场传播至边界处的处理方法采用PML吸收边界条件，震源波场恢复部分采用逆PML波场恢复技术，计算部分采用GPU并行加速技术。本发明采用PML吸收边界条件，比随机边界条件更能够有效的消除边界反射干扰，采用逆PML波场恢复技术，相比于传统的波场存储策略，在不增加计算量的前提下，大幅减少了存储量的需求，而且波场恢复效果良好，无振幅衰减，采用GPU并行加速技术,提高了运算效率。
Description
技术领域
本发明属于逆时偏移领域，尤其涉及一种碳酸盐岩储层的逆时偏移方法。
背景技术
碳酸盐岩背景技术
碳酸盐岩储层结构和渗流通道多以孔洞及将其连接的裂缝构造为主，该类储层通常为高渗透类储层。国内外大量油气勘探实践证明，碳酸盐岩储层主要具备高渗透性，孔隙构造微幅且复杂，非均质性较强等特点，也为地震勘探方法精确识别该类储层带来了较大的困难。由于碳酸盐岩储层构造的反射波散射严重，能量较弱，且构造复杂致使反射波相干性较差，运用传统的偏移方法难以有效识别碳酸盐岩储层。近年来，一些学者致力于碳酸盐岩储层的成像技术研究，张红军(2011)针对盐下碳酸盐岩储层进行了叠前深度偏移成像分析；王小卫(2011)利用波动方程叠前偏移方法对在碳酸盐岩地层成像；白英哲(2011)将Q偏移应用到了塔里木盆地碳酸盐岩储层的成像处理中。但由于碳酸盐岩的缝洞孔隙构造常伴随有由于构造运动所形成的断裂和大角度的形态扭曲，常规的偏移方法难以对其进行精确刻画。
逆时偏移的成像理论基础来自于Claerbout提出的时间一致性成像原理(Claerbout，1971)，即成像点存在于地层内下行波波至时间与某一上行波波至时间相一致之处。对于逆时偏移来说就是成像点位于震源波场外推与接收波场外推时间相一致之处。因此在逆时偏移的算法过程中主要涉及波场的正向及反向延拓计算和激发点波场与检波点波场的相关成像计算。
边界条件方面，现有技术一种是采用随机边界条件，通过在边界外围设置随机减速层，将边界处波场以随机噪音的形式保存下来，因此只需保存最大时刻波场，即可通过方程的逆运算，获取激发点的历史时刻波场，进行逆时成像；另一种方案是针对吸收边界条件，进行波场保存，在逆时成像过程中，通过读取磁盘中记录下的历史时刻波场进行成像计算。针对后者，也可采用设置检查点波场，保存部分历史波场，然后通过调取检查点波场进行正向波场延拓来获取相应时刻的波场信息。
计算效率方面，现有技术可分为串行计算和并行计算两种方向。串行计算方案计算成本较高，影响普及应用。并行计算分为MPI并行和GPU并行加速技术。
现有技术方案在进行有限差分求解波动方程进行波场模拟过程中，采用常规的有限差分系数，会产生一定的数值频散，影响数值模拟的准确性。
通过采用随机边界条件，虽然解决了存储量的问题，但是会产生随机噪音，进一步影响了成像的精度。若采用PML完全匹配层吸收边界条件，则需要对波场信息进行存储，其存储量颇为巨大，是目前计算机技术很难实现的，较为严重的影响了逆时偏移技术的工业应用。而针对PML边界条件结合检查点的存储方案，虽然一定程度上的解决了存储问题，但是对计算量的需求却大大的提升了。
计算效率方面，采用传统的串行计算模式，计算效率比较低。采用MPI并行技术，可以一定程度的提高计算效率，但是需要很大的经济成本。
发明内容
本发明的目的在于提供一种碳酸盐岩储层的逆时偏移方法，旨在解决逆时偏移成本高以及与存储量需求大的问题，同时，有效消除由于数值计算引发的边界噪音，提高成像算法的精度，充分利用逆时偏移算法在复杂介质成像方面的优势，提高碳酸盐岩储层裂缝及溶洞的识别精度。
本发明是这样实现的，一种碳酸盐岩储层的逆时偏移方法采用有限差分方法求解波动方程进行波场模拟，波场传播至边界处的处理方法采用PML吸收边界条件，震源波场恢复部分采用逆PML波场恢复技术，计算部分采用GPU并行加速技术。
进一步，采用有限差分方法求解波动方程进行波场模拟的具体方法为：
如果用U表示某一时刻t二维空间任一点(x，z)处的位移，二维声波方程可表示如下：
式中，U为声波波场，v为速度场，时间导数采用二阶中心差分、空间导数为2N阶差分精度的二维声波方程的高阶有限差分格式为：
式中，
其中，al为有限差分系数，即
进一步，波场传播至边界处的处理方法采用PML吸收边界条件，具体方法为：
在边界外围加入PML吸收衰减层，当波场在有效区域传播时，不发生衰减，当波场进入吸收衰减层后，会发生一定规律的吸收衰减，衰减的形式取决于衰减函数的选取，选择cos衰减函数，加入PML边界条件后的声波方程如下：
其差分格式如下：
cos衰减函数公式如下：
i＝0,1,2···L，式中L为吸收衰减层数，i为计算点位置距吸收层内边界距离，B为衰减常数。
进一步，通过波场的正向模拟后，在震源波场恢复部分，采用逆PML波场恢复技术，其实现方程如下：
对上述的加入PML边界条件后的声波方程进行离散，设定Un-1(x,z)为方程的新数值解，整理差分方程如下：
进一步，所述的碳酸盐岩储层的逆时偏移方法采用GPU并行加速技术，将计算所需数据传至显存内，然后调用编写的GPU版本算法代码，对算法进行加速计算，计算过程中，调用GPU卡内部的高速共享存储器来降低数据的访问延迟，进一步提高计算效率。
效果汇总
本发明采用PML吸收边界条件，比随机边界条件更能够有效的消除边界反射干扰，采用逆PML波场恢复技术，相比于传统的波场存储策略，在不增加计算量的前提下，大幅减少了存储量的需求，而且波场恢复效果良好，无振幅衰减，采用GPU并行加速技术,提高了运算效率。
附图说明
图1是本发明实施例提供的cos函数衰减趋势图；
图2是本发明实施例提供的边界吸收对比效果图；
图3是本发明实施例提供的PML吸收边界逆时偏移结果图；
图4是本发明实施例提供的随机边界逆时偏移结果图；
图5是本发明实施例提供的波场振幅对比图；
图6是本发明实施例提供的逆PML逆时偏移剖面；
图7是本发明实施例提供的checkpoint逆时偏移剖面；
图8是本发明实施例提供的碳酸盐岩构造速度模型；
图9是本发明实施例提供的碳酸盐岩构造偏移剖面。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
本发明是这样实现的，一种碳酸盐岩储层的逆时偏移方法采用有限差分方法求解波动方程进行波场模拟，波场传播至边界处的处理方法采用PML吸收边界条件，震源波场恢复部分采用逆PML波场恢复技术，计算部分采用GPU并行加速技术。
进一步，采用有限差分方法求解波动方程进行波场模拟的具体方法为：
如果用U表示某一时刻t二维空间任一点(x，z)处的位移，二维声波方程可表示如下：
式中，U为声波波场，v为速度场，时间导数采用二阶中心差分、空间导数为2N阶差分精度的二维声波方程的高阶有限差分格式为：
式中，
其中，al为有限差分系数，即
进一步，波场传播至边界处的处理方法采用PML吸收边界条件，具体方法为：
在边界外围加入PML吸收衰减层，当波场在有效区域传播时，不发生衰减，当波场进入吸收衰减层后，会发生一定规律的吸收衰减，衰减的形式取决于衰减函数的选取，选择cos衰减函数，加入PML边界条件后的声波方程如下：
其差分格式如下：
cos衰减函数公式如下：
i＝0,1,2···L，式中L为吸收衰减层数，i为计算点位置距吸收层内边界距离，B为衰减常数。
图1为衰减函数A的趋势曲线，通过这种衰减趋势的设置，可以使得波场被逐步递进的吸收衰减，可以有效避免波场因突然变化而产生强反射。
图2为吸收效果对比，左上部分加入了PML吸收边界，波场传播至吸收层后，被逐步吸收，右下部分未加入边界条件，波场发生反射干扰。
通过图3和图4的对比可以看出，采用随机边界条件进行逆时偏移，会在浅层部分产生一定的随机噪音，而采用PML吸收边界条件，随机噪音会得以的消除。
进一步，通过波场的正向模拟后，在震源波场恢复部分，采用逆PML波场恢复技术，其实现方程如下：
对上述的加入PML边界条件后的声波方程进行离散，设定Un-1(x,z)为方程的新数值解，整理差分方程如下：
由于该方程的不具备持续稳定性，所以在保证逆PML方程稳定性的前提下，保存尽量少的波场信息，进而实现波场恢复。下图是模拟波场与逆PML波场恢复的振幅对比图，从图4中黑色和白色点位的重合效果可以看出，恢复的波场振幅与原始波场振幅完全吻合，没有任何振幅损失。
图6所示的是采用本发明的逆PML波场恢复技术进行逆时偏移的成像剖面，图7所示的是checkpoint波场存储进行叠前逆时偏移的成像剖面，通过对比可以看出二者均能够有效的压制边界噪音，且二者的成像质量无明显差异。
进一步，由于该算法的计算量的需求较为巨大，因此需要对计算效率进行提升，采用GPU并行加速技术，将计算所需数据传至显存内，然后调用编写的GPU版本算法代码，对算法进行加速计算，计算过程中，调用GPU卡内部的高速共享存储器来降低数据的访问延迟，进一步提高计算效率。
图8和图9为应用本发明方法进行的逆时偏移试算，结果表明，应用本发明方案可以有效的对碳酸盐岩类精细构造进行成像，效果良好。
本发明采用PML吸收边界条件，比随机边界条件更能够有效的消除边界反射干扰，采用逆PML波场恢复技术，相比于传统的波场存储策略，在不增加计算量的前提下，大幅减少了存储量的需求，而且波场恢复效果良好，无振幅衰减，采用GPU并行加速技术,提高了运算效率。
上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性的劳动即可做出的各种修改或变形仍在本发明的保护范围之内。
Claims (1)
1.一种碳酸盐岩储层的逆时偏移方法，其特征在于，所述的碳酸盐岩储层的逆时偏移方法采用有限差分方法求解波动方程进行波场模拟，波场传播至边界处的处理方法采用PML吸收边界条件，震源波场恢复部分采用逆PML波场恢复技术，计算部分采用GPU并行加速技术；
该方法采用有限差分方法求解波动方程进行波场模拟，具体方法为：
用U表示某一时刻t二维空间任一点(x，z)处的位移，二维声波方程可表示如下：
<mrow>
<mfrac>
<mrow>
<msup>
<mo>&part;</mo>
<mn>2</mn>
</msup>
<mi>U</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<msup>
<mi>x</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>+</mo>
<mfrac>
<mrow>
<msup>
<mo>&part;</mo>
<mn>2</mn>
</msup>
<mi>U</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<msup>
<mi>z</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<msup>
<mi>v</mi>
<mn>2</mn>
</msup>
</mfrac>
<mfrac>
<mrow>
<msup>
<mo>&part;</mo>
<mn>2</mn>
</msup>
<mi>U</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<msup>
<mi>t</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
</mrow>
式中，U为声波波场，v为速度场，时间导数采用二阶中心差分、空间导数为2N阶差分精度的二维声波方程的高阶有限差分格式为：
<mrow>
<msup>
<mi>U</mi>
<mrow>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mn>2</mn>
<msup>
<mi>U</mi>
<mi>n</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msup>
<mi>U</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msup>
<mi>dt</mi>
<mn>2</mn>
</msup>
<msup>
<mi>v</mi>
<mn>2</mn>
</msup>
<mo>{</mo>
<msubsup>
<mi>C</mi>
<mi>x</mi>
<mn>2</mn>
</msubsup>
<mo>&lsqb;</mo>
<msup>
<mi>U</mi>
<mi>n</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>+</mo>
<msubsup>
<mi>C</mi>
<mi>z</mi>
<mn>2</mn>
</msubsup>
<mo>&lsqb;</mo>
<msup>
<mi>U</mi>
<mi>n</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>}</mo>
</mrow>
式中，
<mrow>
<msubsup>
<mi>C</mi>
<mi>x</mi>
<mn>2</mn>
</msubsup>
<mo>&lsqb;</mo>
<msup>
<mi>U</mi>
<mi>n</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>l</mi>
<mo>=</mo>
<mo>-</mo>
<mi>N</mi>
</mrow>
<mi>N</mi>
</munderover>
<mfrac>
<msub>
<mi>a</mi>
<mi>l</mi>
</msub>
<mrow>
<msup>
<mi>&Delta;x</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mi>U</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>+</mo>
<mi>l</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msubsup>
<mi>C</mi>
<mi>z</mi>
<mn>2</mn>
</msubsup>
<mo>&lsqb;</mo>
<msup>
<mi>U</mi>
<mi>n</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>l</mi>
<mo>=</mo>
<mo>-</mo>
<mi>N</mi>
</mrow>
<mi>N</mi>
</munderover>
<mfrac>
<msub>
<mi>a</mi>
<mi>l</mi>
</msub>
<mrow>
<msup>
<mi>&Delta;z</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mi>U</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>+</mo>
<mi>l</mi>
<mo>)</mo>
</mrow>
</mrow>
其中，al为有限差分系数，即
<mrow>
<msub>
<mi>a</mi>
<mi>l</mi>
</msub>
<mo>=</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mo>)</mo>
</mrow>
<mrow>
<mi>l</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<msup>
<mi>l</mi>
<mn>2</mn>
</msup>
</mfrac>
<mfrac>
<mrow>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mo>)</mo>
</mrow>
<mrow>
<mi>l</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<munderover>
<mi>&Pi;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mi>i</mi>
<mo>&NotEqual;</mo>
<mi>l</mi>
</mrow>
<mi>N</mi>
</munderover>
<msup>
<mi>i</mi>
<mn>2</mn>
</msup>
</mrow>
<mrow>
<munderover>
<mi>&Pi;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>l</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</munderover>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>l</mi>
<mn>2</mn>
</msup>
<mo>-</mo>
<msup>
<mi>i</mi>
<mn>2</mn>
</msup>
</mrow>
<mo>)</mo>
</mrow>
<munderover>
<mi>&Pi;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mi>l</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>i</mi>
<mn>2</mn>
</msup>
<mo>-</mo>
<msup>
<mi>l</mi>
<mn>2</mn>
</msup>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>;</mo>
</mrow>
其中，波场传播至边界处的处理方法采用PML吸收边界条件，具体方法为：
在边界外围加入PML吸收衰减层，当波场在边界内区域传播时，不发生衰减，当波场进入吸收衰减层后，会发生一定规律的吸收衰减，衰减的形式取决于衰减函数的选取，选择cos衰减函数，加入PML边界条件后的声波方程如下：
<mrow>
<mfrac>
<mrow>
<msup>
<mo>&part;</mo>
<mn>2</mn>
</msup>
<mi>U</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<msup>
<mi>t</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>+</mo>
<mi>A</mi>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>&CenterDot;</mo>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>U</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>t</mi>
</mrow>
</mfrac>
<mo>+</mo>
<msup>
<mi>A</mi>
<mn>2</mn>
</msup>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>&CenterDot;</mo>
<mi>U</mi>
<mo>=</mo>
<msup>
<mi>v</mi>
<mn>2</mn>
</msup>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<msup>
<mo>&part;</mo>
<mn>2</mn>
</msup>
<mi>U</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<msup>
<mi>x</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>+</mo>
<mfrac>
<mrow>
<msup>
<mo>&part;</mo>
<mn>2</mn>
</msup>
<mi>U</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<msup>
<mi>z</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
</mrow>
其差分格式如下：
<mrow>
<msup>
<mi>U</mi>
<mrow>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mrow>
<mn>2</mn>
<mo>-</mo>
<msup>
<mi>dt</mi>
<mn>2</mn>
</msup>
<msup>
<mi>A</mi>
<mn>2</mn>
</msup>
</mrow>
<mrow>
<mn>1</mn>
<mo>+</mo>
<mi>t</mi>
<mi>A</mi>
</mrow>
</mfrac>
<msup>
<mi>U</mi>
<mi>n</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<mrow>
<mn>1</mn>
<mo>-</mo>
<mi>t</mi>
<mi>A</mi>
</mrow>
<mrow>
<mn>1</mn>
<mo>+</mo>
<mi>t</mi>
<mi>A</mi>
</mrow>
</mfrac>
<msup>
<mi>U</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<mrow>
<msup>
<mi>dt</mi>
<mn>2</mn>
</msup>
<msup>
<mi>v</mi>
<mn>2</mn>
</msup>
</mrow>
<mrow>
<mn>1</mn>
<mo>+</mo>
<mi>t</mi>
<mi>A</mi>
</mrow>
</mfrac>
<mo>{</mo>
<msubsup>
<mi>C</mi>
<mi>x</mi>
<mn>2</mn>
</msubsup>
<mo>&lsqb;</mo>
<msup>
<mi>U</mi>
<mi>n</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>+</mo>
<msubsup>
<mi>C</mi>
<mi>z</mi>
<mn>2</mn>
</msubsup>
<mo>&lsqb;</mo>
<msup>
<mi>U</mi>
<mi>n</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>}</mo>
</mrow>
cos衰减函数公式如下：
<mrow>
<mi>A</mi>
<mo>=</mo>
<mi>B</mi>
<mo>{</mo>
<mn>1</mn>
<mo>-</mo>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mo>&lsqb;</mo>
<mfrac>
<mrow>
<mi>&pi;</mi>
<mrow>
<mo>(</mo>
<mi>l</mi>
<mo>-</mo>
<mi>i</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mn>2</mn>
<mi>L</mi>
</mrow>
</mfrac>
<mo>&rsqb;</mo>
<mo>}</mo>
<mo>,</mo>
<mi>i</mi>
<mo>=</mo>
<mn>0</mn>
<mo>,</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>...</mo>
<mi>L</mi>
<mo>,</mo>
</mrow>
式中L为吸收衰减层数，i为计算点位置距吸收层内边界距离，B为衰减常数；
其中，通过波场的正向模拟后，在震源波场恢复部分，采用逆PML波场恢复技术，实现方程如下：
加入PML边界条件后的声波方程
<mrow>
<mfrac>
<mrow>
<msup>
<mo>&part;</mo>
<mn>2</mn>
</msup>
<mi>U</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<msup>
<mi>t</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>-</mo>
<mi>A</mi>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>&CenterDot;</mo>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>U</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>t</mi>
</mrow>
</mfrac>
<mo>-</mo>
<msup>
<mi>A</mi>
<mn>2</mn>
</msup>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>&CenterDot;</mo>
<mi>U</mi>
<mo>=</mo>
<msup>
<mi>v</mi>
<mn>2</mn>
</msup>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<msup>
<mo>&part;</mo>
<mn>2</mn>
</msup>
<mi>U</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<msup>
<mi>x</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>+</mo>
<mfrac>
<mrow>
<msup>
<mo>&part;</mo>
<mn>2</mn>
</msup>
<mi>U</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<msup>
<mi>z</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
</mrow>
进行离散，设定Un-1(x,z)为方程的新数值解，整理差分方程如下：
<mrow>
<msup>
<mi>U</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mrow>
<mn>2</mn>
<mo>-</mo>
<msup>
<mi>dt</mi>
<mn>2</mn>
</msup>
<msup>
<mi>A</mi>
<mn>2</mn>
</msup>
</mrow>
<mrow>
<mn>1</mn>
<mo>-</mo>
<mi>t</mi>
<mi>A</mi>
</mrow>
</mfrac>
<msup>
<mi>U</mi>
<mi>n</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<mrow>
<mn>1</mn>
<mo>+</mo>
<mi>t</mi>
<mi>A</mi>
</mrow>
<mrow>
<mn>1</mn>
<mo>-</mo>
<mi>t</mi>
<mi>A</mi>
</mrow>
</mfrac>
<msup>
<mi>U</mi>
<mrow>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<mrow>
<msup>
<mi>dt</mi>
<mn>2</mn>
</msup>
<msup>
<mi>v</mi>
<mn>2</mn>
</msup>
</mrow>
<mrow>
<mn>1</mn>
<mo>-</mo>
<mi>t</mi>
<mi>A</mi>
</mrow>
</mfrac>
<mo>{</mo>
<msubsup>
<mi>C</mi>
<mi>x</mi>
<mn>2</mn>
</msubsup>
<mo>&lsqb;</mo>
<msup>
<mi>U</mi>
<mi>n</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>+</mo>
<msubsup>
<mi>C</mi>
<mi>z</mi>
<mn>2</mn>
</msubsup>
<mo>&lsqb;</mo>
<msup>
<mi>U</mi>
<mi>n</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>}</mo>
<mo>;</mo>
</mrow>
其中，所述的碳酸盐岩储层的逆时偏移方法采用GPU并行加速技术，将计算所需数据传至显存内，然后调用编写的GPU版本算法代码，对算法进行加速计算，计算过程中，调用GPU卡内部的高速共享存储器来降低数据的访问延迟，进一步提高计算效率。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410325458.4A CN104133987B (zh) | 2014-07-09 | 2014-07-09 | 一种碳酸盐岩储层的逆时偏移方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410325458.4A CN104133987B (zh) | 2014-07-09 | 2014-07-09 | 一种碳酸盐岩储层的逆时偏移方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104133987A CN104133987A (zh) | 2014-11-05 |
| CN104133987B true CN104133987B (zh) | 2018-01-09 |
Family
ID=51806662
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410325458.4A Expired - Fee Related CN104133987B (zh) | 2014-07-09 | 2014-07-09 | 一种碳酸盐岩储层的逆时偏移方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104133987B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106354897A (zh) * | 2015-07-17 | 2017-01-25 | 中国石油化工股份有限公司 | 一种基于gpu的卷积最佳匹配层边界条件实现方法 |
| CN105403919B (zh) * | 2015-11-11 | 2018-02-02 | 中国石油天然气集团公司 | 一种逆时偏移成像方法及装置 |
| CN107102353B (zh) * | 2017-05-08 | 2019-09-03 | 厦门大学 | 基于高阶差分方法的弹性波方程逆时偏移成像方法 |
| CN108051855B (zh) * | 2017-12-13 | 2019-02-15 | 国家深海基地管理中心 | 一种基于拟空间域声波方程的有限差分计算方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103142327A (zh) * | 2013-03-20 | 2013-06-12 | 张蕾 | 实验鼠无创固定器 |
| CN103777238A (zh) * | 2012-10-17 | 2014-05-07 | 中国石油化工股份有限公司 | 一种纯纵波各向异性波场模拟方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8614930B2 (en) * | 2011-03-23 | 2013-12-24 | Chevron U.S.A. Inc. | System and method for seismic data modeling and migration |
-
2014
- 2014-07-09 CN CN201410325458.4A patent/CN104133987B/zh not_active Expired - Fee Related
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103777238A (zh) * | 2012-10-17 | 2014-05-07 | 中国石油化工股份有限公司 | 一种纯纵波各向异性波场模拟方法 |
| CN103142327A (zh) * | 2013-03-20 | 2013-06-12 | 张蕾 | 实验鼠无创固定器 |
Non-Patent Citations (4)
| Title |
|---|
| 地震叠前逆时偏移的有效边界存储策略;王保利 等;《地球物理学报》;20120715;第55卷(第7期);2412-2421页 * |
| 基于GPU并行加速的叠前逆时偏移方法;石颖 等;《东北石油大学学报》;20120926;第36卷(第4期);111-116页 * |
| 完全匹配层吸收边界条件应用研究;王维红 等;《地球物理学进展》;20131015;第28卷(第5期);2508-2015页 * |
| 碳酸盐岩孔缝洞复杂模型逆时偏移成像方法研究;黄建平 等;《中国地球物理2011》;20111017;578页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104133987A (zh) | 2014-11-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104133987B (zh) | 一种碳酸盐岩储层的逆时偏移方法 | |
| CN103424777B (zh) | 一种提高地震成像分辨率的方法 | |
| CN102269820B (zh) | 一种三维地震叠前逆时偏移成像方法 | |
| CN102590859B (zh) | 垂向各向异性介质准p波方程逆时偏移方法 | |
| CN108181653B (zh) | 针对vti介质逆时偏移方法、设备及介质 | |
| CN106526674A (zh) | 一种三维全波形反演能量加权梯度预处理方法 | |
| CN105549087B (zh) | 一种煤矿井下槽波地震勘探的走时和振幅联合反演方法 | |
| CN106526677A (zh) | 一种海上自适应压制鬼波的宽频逆时偏移成像方法 | |
| CN104991269A (zh) | 一种边缘引导和结构约束的全波形反演快速方法 | |
| CN102116869A (zh) | 高精度叠前域最小二乘偏移地震成像技术 | |
| CN103995288A (zh) | 一种高斯束叠前深度偏移方法及装置 | |
| CN103675915B (zh) | 基于地震资料估计地层横向相对品质因子的方法和装置 | |
| CN103926623A (zh) | 一种压制逆时偏移低频噪音的方法 | |
| CN103489159A (zh) | 基于三边结构导向滤波的三维地震数据图像降噪方法 | |
| CN109001813A (zh) | 一种压制多次波的方法、装置及系统 | |
| CN108051855A (zh) | 一种基于拟空间域声波方程的有限差分计算方法 | |
| CN106249290A (zh) | 一种利用多级数据融合建立表层速度结构模型的方法 | |
| CN102565852B (zh) | 针对储层含油气性检测的角度域叠前偏移数据处理方法 | |
| CN105447225A (zh) | 一种应用于声波有限差分数值模拟的组合吸收边界条件 | |
| CN107918153A (zh) | 一种地震信号相干性高精度检测方法 | |
| Liu et al. | The series solution to the modified mild-slope equation for wave scattering by Homma islands | |
| CN105353409B (zh) | 一种用于抑制全波形反演震源编码串扰噪音的方法和系统 | |
| CN104280768B (zh) | 一种适用于逆时偏移的吸收边界条件方法 | |
| CN106199692A (zh) | 一种基于gpu的波动方程反偏移方法 | |
| Zhou et al. | An iterative factored topography-dependent eikonal solver for anisotropic media |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20180109 Termination date: 20190709 |