CN108181653B - 针对vti介质逆时偏移方法、设备及介质 - Google Patents
针对vti介质逆时偏移方法、设备及介质 Download PDFInfo
- Publication number
- CN108181653B CN108181653B CN201810039657.7A CN201810039657A CN108181653B CN 108181653 B CN108181653 B CN 108181653B CN 201810039657 A CN201810039657 A CN 201810039657A CN 108181653 B CN108181653 B CN 108181653B
- Authority
- CN
- China
- Prior art keywords
- wave
- continuation
- reverse
- wave field
- equation
- 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
- 238000000034 method Methods 0.000 title claims abstract description 75
- 238000013508 migration Methods 0.000 title claims abstract description 75
- 230000005012 migration Effects 0.000 title claims abstract description 75
- 238000003384 imaging method Methods 0.000 claims abstract description 37
- 238000004364 calculation method Methods 0.000 claims abstract description 13
- 238000003860 storage Methods 0.000 claims description 24
- 238000004088 simulation Methods 0.000 claims description 22
- 238000004590 computer program Methods 0.000 claims description 12
- 239000000284 extract Substances 0.000 claims description 6
- 239000002360 explosive Substances 0.000 claims description 4
- 238000010521 absorption reaction Methods 0.000 claims description 3
- 239000003989 dielectric material Substances 0.000 abstract description 7
- 238000012545 processing Methods 0.000 description 11
- 238000010586 diagram Methods 0.000 description 10
- 230000008569 process Effects 0.000 description 7
- 239000006185 dispersion Substances 0.000 description 6
- 238000004891 communication Methods 0.000 description 5
- 230000001154 acute effect Effects 0.000 description 4
- 238000004422 calculation algorithm Methods 0.000 description 4
- 238000013499 data model Methods 0.000 description 4
- 241000208340 Araliaceae Species 0.000 description 3
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 3
- 235000003140 Panax quinquefolius Nutrition 0.000 description 3
- 229910003460 diamond Inorganic materials 0.000 description 3
- 239000010432 diamond Substances 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000013213 extrapolation Methods 0.000 description 3
- 235000008434 ginseng Nutrition 0.000 description 3
- 238000004519 manufacturing process Methods 0.000 description 3
- 230000000644 propagated effect Effects 0.000 description 3
- 230000005856 abnormality Effects 0.000 description 2
- 230000001133 acceleration Effects 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 239000003208 petroleum Substances 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000000151 deposition Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 239000000835 fiber Substances 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000011514 reflex Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
- 238000010200 validation analysis Methods 0.000 description 1
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)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供的针对VTI介质逆时偏移方法、设备及介质，方法为：在炮集数据中提取震源坐标，并生成震源子波；根据正演方案，通过qP波波动方程对震源子波沿时间正向延拓，在延拓过程中选择若干部分波场进行存储；通过qP波波动方程利用若干部分波场沿时间逆向延拓，得到震源波场；在炮集数据中提取地震数据，通过qP波波动方程对插值后的地震数据沿时间逆向延拓，得到检波点波场；将震源波场和检波点波场进行互相关成像，得到一炮偏移剖面，重复上述步骤，得到多炮偏移剖面；将多个炮偏移剖面进行叠加，得到最终成像剖面。本发明在大幅减小计算量的情况下，能够获得相同运动学特征的纵波分量，使成像剖面的反射同相轴准确归位，提升复杂介质的成像精度。
Description
技术领域
本发明涉及地震勘测技术领域，尤其涉及针对VTI介质逆时偏移方法及设备领域。
背景技术
地震勘探是寻找油气资源的主要手段之一，在勘探地震学研究领域，早期由Whitmore(1983)提出的逆时偏移方法是获取地下构造形态最有效的手段之一。近年来，由于计算机科学技术的长足进步，为逆时偏移方法带来了工业化的可能，利用有限差分数值解法实现了基于声波方程各向同性介质的逆时偏移方法，油气工业界将该方法应用到实际生产中，为油气的增产上储提供有力的技术支撑。随着油气勘探开发的不断深入，可供开采传统的一些构造简单的油气藏逐渐减少，因此非常规油气资源在能源格局中的地位愈发重要，而非常规油气储层多呈现出具有垂直对称轴的横向各向同性特征，即地震波的传播速度随传播角度的变化而变化。若针对VTI介质仍采用基于各向同性介质偏移理论，将出现地震波的振幅和相位不准确等问题，导致复杂介质的成像精度难以保证。因此，针对VTI介质逆时偏移方法在油气勘探中彰显了更大的研究潜力和发展前景。
在理论上，虽然对于VTI介质基于弹性波方程的逆时偏移理论是最佳的方案，并且对于模型数据可以获得理想的成像剖面，但是在实施过程中会遇到诸多问题，如纵横波耦合的波场会为成像剖面带来串扰噪音，而弹性波波场分离增加了逆时偏移的难度；如地下介质的弹性参数矩阵在实际生产中难以获取，为弹性波波场的数值模拟带来了未知的因素；再如求解弹性波过程中，由于横波的速度小于纵波速度，考虑波场传播的稳定性因素，在数值模拟过程中只能采用较小的空间网格，这无形地增加了计算量。所以基于弹性波方程的逆时偏移理论并不适用于地震勘探的需求。
时至今日，地震勘探仍然以纵波勘探为主，因此为了将VTI介质逆时偏移方法从理论走向工业应用，高效地实现准纵(qP)波在VTI介质中精确延拓成为该方法的核心问题。通过声学近似方法提出适用于VTI介质四阶的波动方程，开启了针对该介质qP波场精确描述的新篇章。通过引入辅助变量将四阶波动方程转换为两个耦合的波动方程，随后地球科学工作者在四阶波动方程中引入不同的辅助变量得到一些列二阶耦合波动方程，这极大地简化了数值模拟复杂度，但是在各向异性参数变化剧烈的情况下，利用现有的这些波动方程对波场进行延拓将导致波场异常，而无法应用到逆时偏移方法中。
在逆时偏移过程中，若采用互相关成像条件，同时需获取每个时刻的震源波场和检波点波场。利用反延拓算子将地震记录加载到检波点直接求取检波点波场，对于获得震源波场主要有两种途径，一种是在数值模拟震源波场的过程中，将每个时刻的震源波场保存至计算机内存中，互相关成像时从计算机内存中读取，虽然这种方式简单易行，但是将占用大量的计算机内存；另外一种方式是通过保存部分震源波场信息实现震源波场重建，该种方式能够有效地减小内存的占用，并且得到了理想的重建效果。而对于大数据模型的VTI介质而言第一种方式无法满足内存占用过大的问题，即使利用第二种方式也另计算硬件掣肘。
在基于有限差分数值解法的逆时偏移计算中，波场延拓占用主要的机时。对于小数据模型的逆时偏移计算，考虑到计算计算技术的准入门槛、计算设备的成本和计算耗时等因素，图形处理器设备是一个理想的选择。而对于复杂的大数据模型的逆时偏移计算，由于中央处理器设备在架构设计上，与算数单元相比逻辑控制单元占主体部分，更注重逻辑关系的判断，所以中央处理器设备的计算能力已经严重影响数据的处理能力，计算周期甚至是难以令人接受的，利用中央处理器设备实现逆时偏移算法的计算耗时成为其投入工业应用的瓶颈。
因此，现有技术中，利用全弹性波方程进行偏移计算量大而无法满足工业化生产实际需求；在各向异性参数变化剧烈的情况下，利用声学近似方法得到的准纵波方程求解波场时存在不稳定现象；针对大模型数据时，现有的存储策略无法满足图形处理器对存储的要求。
发明内容
针对上述技术问题，本发明提供一种针对VTI介质逆时偏移方法及设备，利用准纵波方程实现地震波在VTI介质中的数值模拟，与利用弹性波方程相比准纵波方程在大幅减小计算量的情况下，能够获得相同运动学特征的纵波分量，基于qP波方程建立针对VTI介质逆时偏移成像的处理流程，使成像剖面的反射同相轴准确归位，提升复杂介质的成像精度。
为解决上述技术问题，本发明提供的技术方案是：
第一方面，本发明提供一种针对VTI介质逆时偏移方法，包括：
步骤S1，获取炮集数据；
步骤S2，在所述炮集数据中提取震源坐标，并生成震源子波；
步骤S3，根据正演方案，将震源子波加载至波场中，得到震源波场，通过qP波波动方程对所述震源波场沿时间正向延拓，共进行Nc次延拓，Nc为t0至tNc时间段内设定的延拓次数，并在延拓过程中选择若干部分波场进行存储；
步骤S4，根据所述正演方案，通过所述qP波波动方程对所述若干部分波场沿时间分段逆向延拓，共进行Nc次延拓，实现对所述震源波场的重建；
步骤S5，在所述炮集数据中提取地震数据，并对所述地震数据进行插值，得到插值后的地震数据；
步骤S6，根据所述正演方案，通过所述qP波波动方程对所述插值后的地震数据沿时间逆向延拓，共进行Nc次延拓，得到检波点波场；
步骤S7，将所述震源波场和所述检波点波场进行互相关成像，得到一炮偏移剖面，重复上述步骤，得到多炮偏移剖面；
步骤S8，将所述多炮偏移剖面进行叠加，得到最终成像剖面。
本发明提供的针对VTI介质逆时偏移方法，利用准纵波方程实现地震波在VTI介质中的数值模拟，与利用弹性波方程相比准纵波方程在大幅减小计算量的情况下，能够获得相同运动学特征的纵波分量，基于qP波方程建立针对VTI介质逆时偏移成像的处理流程，使成像剖面的反射同相轴准确归位，提升复杂介质的成像精度。
进一步地，所述步骤S2之后，还包括：
对所述震源子波进行插值，得到插值后的震源子波。
进一步地，基于qP波波动方程进行数值模拟计算，得到模拟地震波场，具体为：
利用高阶有限差分数值方法求解qP波波动方程，采用PML完全匹配层边界条件吸收边界反射，得到模拟地震波场。
进一步地，在所述qP波波动方程中引入稳定性因子σ；
利用稳定性因子公式实现对横波沿着对称轴方向的传播速度Vsz的计算，其中，所述稳定性因子σ为0.75，ε和δ为Thomsen参数，Vpz为纵波沿着对称轴方向的传播速度。
进一步地，所述步骤S3中，根据正演方案，通过qP波波动方程对所述震源子波沿时间正向延拓，具体为：
根据所述纵波和横波沿着对称轴方向的传播速度以及Thomsen参数作为基于qP波波动方程的VTI介质逆时偏移成像方法的模型参数；
将地震子波作为激发震源生成地震波场，利用有限差分数值解法，对所述震源波场沿时间正向Nc次延拓，Nc为t0至tNc时间段内设定的延拓次数；
其中，qP波波动方程如以下公式所示，
式中，Vpz和Vsz分别为纵波和横波沿着对称轴方向的传播速度，ε和δ为Thomsen参数，Vpx为纵波沿着垂直于对称轴方向的传播速度，Vpn为动校正速度，p为声压，q为辅助变量。
进一步地，所述通过基于有效边界的检查点存储策略进行若干部分波场的存储，具体为：
假设所述震源波场正向延拓的总时间为tNc，在所述震源波场中确定N个检查点；
在所述震源波场正向延拓的过程中，保存每个检查点时刻及前一时刻的p和q全部波场；
将所述每个检查点时刻及前一时刻的p和q全部波场作为若干部分波场进行存储，其中，p为声压，q为辅助变量。
进一步地，所述步骤S4，具体为：
在对所述检查点N-1至检查点N时段内波场进行重建时，读取所述检查点N-1时刻及前一时刻的p和q全部波场信息作为初值进行c次正向延拓，每次延拓记录内边界层信息，保存至图形处理器的显存p和q缓冲区，所述显存p和q缓冲区为两个检查点之间p和q内边界大小；
将tNc和tNc-1两时刻的p和q全部波场作为反向延拓的初值，每延拓一次，将所述保存至显存上p和q缓冲区对应时刻内边界信息载入重建波场中，反向延拓c-1次得到重建t(N-1)c时刻波场；
重复上述步骤重建tNc至t0波场，得到震源波场。
进一步地，在进行所述正向延拓和逆向延拓时，通过图形处理器进行计算。
第二方面，本发明提供一种针对VTI介质逆时偏移设备，包括：至少一个处理器、至少一个存储器以及存储在所述存储器中的计算机程序指令，当所述计算机程序指令被所述处理器执行时实现如第一方面所述的方法。
第三方面，本发明提供一种计算机可读存储介质，其上存储有计算机程序指令，当所述计算机程序指令被处理器执行时实现如如第一方面所述的方法。
与现有技术相比，本发明的优点：
1、常规声波方程只能对地震波在各向同性介质进行精确描述，本发明利用准纵波方程实现地震波在VTI介质中的数值模拟，与利用弹性波相比准纵波方程在大幅减小计算量的情况下，能够获得相同运动学特征的纵波分量。
2、在准纵波方程中加入稳定性因子，该因子能够能够抑制波场频散，确保地震波在剧烈变化的VTI介质中能够稳定地传播，为获取高精度逆时偏移成像剖面提供前提保障。
3、对于波场的有限差分数值解法，从时间和空间两个角度建立一种检查点存储策略，从而能够极大地节约波场存储空间，同时该策略能降低传统波场重建带来的误差积累，降低逆时偏移算法对计算机硬件的依赖，为方法的工业化带来更广阔的空间。
4、将图形处理器(GPU)设备应用到波场延拓计算中，利用中央处理器处理器(CPU)设备逻辑判断能力强和GPU设备的计算效率高的特点，形成CPU/GPU协同加速技术，从而大幅度减小计算周期。
5、基于qP波方程建立针对VTI介质逆时偏移成像的处理流程，使成像剖面的反射同相轴准确归位，提升复杂介质的成像精度。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍。
图1示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法的流程图；
图2示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中波场存储示意图；
图3a示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中检查点第一存储示意图；
图3b示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中检查点第二存储示意图；
图4a示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中300ms的第一波场图像；
图4b示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中300ms的第二波场图像；
图4c示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中300ms的第三波场图像；
图4d示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中600ms的波场图像；
图5示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中速度模型的示意图；
图6a示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中第一波场重建结果示意图；
图6b示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中第二波场重建结果示意图；
图6c示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中波场重建结果处理示意图；
图7a示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中第一复杂介质逆时偏移成像示意图；
图7b示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中第二复杂介质逆时偏移成像示意图；
图7c示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中第三复杂介质逆时偏移成像示意图；
图7d示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法中逆时偏移成像剖面示意图；
图8示出了本发明实施例所提供的一种针对VTI介质逆时偏移设备的硬件结构示意图。
具体实施方式
下面将结合附图对本发明技术方案的实施例进行详细的描述。以下实施例仅用于更加清楚地说明本发明的技术方案，因此只是作为示例，而不能以此来限制本发明的保护范围。
实施例
图1示出了本发明实施例所提供的一种针对VTI介质逆时偏移方法的流程图；如图1所示，本实施例提供的一种针对VTI介质逆时偏移方法，包括：
步骤S1，获取炮集数据；
步骤S2，在所述炮集数据中提取震源坐标，并生成震源子波；
步骤S3，根据正演方案，将震源子波加载至波场中，得到震源波场，通过qP波波动方程对所述震源波场沿时间正向延拓，共进行Nc次延拓，Nc为t0至tNc时间段内设定的延拓次数，并在延拓过程中选择若干部分波场进行存储；
步骤S4，根据所述正演方案，通过所述qP波波动方程对所述若干部分波场沿时间分段逆向延拓，共进行Nc次延拓，实现对所述震源波场的重建；
步骤S5，在所述炮集数据中提取地震数据，并对所述地震数据进行插值，得到插值后的地震数据；
步骤S6，根据所述正演方案，通过所述qP波波动方程对所述插值后的地震数据沿时间逆向延拓，共进行Nc次延拓，得到检波点波场；
步骤S7，将所述震源波场和所述检波点波场进行互相关成像，得到一炮偏移剖面，重复上述步骤，得到多炮偏移剖面；
步骤S8，将所述多炮偏移剖面进行叠加，得到最终成像剖面。
本发明提供的针对VTI介质逆时偏移方法，具体为针对具有垂直对称轴的横向各向同性(Transversely Isotropy with A Vertical Axis of Symmetry，VTI)介质逆时偏移方法。利用准纵波方程实现地震波在VTI介质中的数值模拟，与利用弹性波相比准纵波方程在大幅减小计算量的情况下，能够获得相同运动学特征的纵波分量，基于qP波方程建立针对VTI介质逆时偏移成像的处理流程，使成像剖面的反射同相轴准确归位，提升复杂介质的成像精度。
优选地，所述步骤S2之后，还包括：
对所述震源子波进行插值，得到插值后的震源子波。
炮集数据中包含不规则的数据，为了使数据规则，对所述震源子波进行插值处理。
本发明的技术方案主要包括两部分，第一部分是正演方案，第二部分是利用qP波波动方程对波场进行延拓，实现VTI介质逆时偏移的关键在于qP波场正演模拟计算，该方程如式1所示，这里采用有限差分法对其求解。
式中，Vpz和Vsz分别为纵波和横波沿着对称轴方向的传播速度，ε和δ为Thomsen参数，Vpx为纵波沿垂直对称轴方向的速度，Vpn为动校正速度，p为声压，q为辅助变量。
其中，采用时间二阶和空间十二阶进行有限差分数值计算。通过上述正演方案进行正向延拓和反向延拓。
优选地，基于qP波波动方程进行数值模拟计算，得到模拟地震波场，具体为：
利用高阶有限差分数值方法求解qP波波动方程，采用PML完全匹配层边界条件吸收边界反射，得到模拟地震波场。
优选地，在所述qP波波动方程中引入稳定性因子σ；
利用稳定性因子公式实现对横波沿着对称轴方向的传播速度Vsz的计算，其中，所述稳定性因子σ为0.75，ε和δ为Thomsen参数，Vpz为纵波沿着对称轴方向的传播速度。
将公式(1)中的Vsz设置为0，即横波沿着对称轴方向的速度为0，在一定程度上可以减少计算量，但是仍然存在沿着非对称轴方向低振幅、低波速的伪横波，这些伪横波的波形特征为菱形交叉状，当地下介质的各向异性参数变化剧烈时，这些交叉状的伪横波将产生严重的频散导致波场异常，从而无法实现逆时偏移算法。为此本发明引入稳定性因子σ，经数值实现验证，σ参数在很大程度上决定qSV波的运动学特征。当设置σ＝0.75时，频散得以压制，确保波场的稳定传播。
在准纵波方程中加入横波速度，该因子能够消除数值模拟中波场出现的频散现象，确保地震波在剧烈变化的VTI介质中能够稳定地延拓，为获取高精度逆时偏移成像剖面提供前提保障。
优选地，所述步骤S3中，根据正演方案，通过qP波波动方程对所述震源子波沿时间正向延拓，具体为：
根据所述纵波和横波沿着对称轴方向的传播速度以及Thomsen参数作为基于qP波波动方程的VTI介质逆时偏移成像方法的模型参数；
将地震子波作为激发震源生成地震波场，利用有限差分数值解法，对所述震源波场沿时间正向Nc次延拓，Nc为t0至tNc时间段内设定的延拓次数；
其中，qP波波动方程如以下公式所示，
式中，Vpz和Vsz分别为纵波和横波沿着对称轴方向的传播速度，ε和δ为Thomsen参数，Vpx为纵波沿着垂直于对称轴方向的传播速度，Vpn为动校正速度，p为声压，q为辅助变量。
优选地，所述通过基于有效边界的检查点存储策略进行若干部分波场的存储，具体为：
假设所述震源波场正向延拓的总时间为tNc，在所述震源波场中确定N个检查点；
在所述震源波场正向延拓的过程中，保存每个检查点时刻及前一时刻的p和q全部波场；
将所述每个检查点时刻及前一时刻的p和q全部波场作为若干部分波场进行存储，其中，p为声压，q为辅助变量。
在所述震源波场正传的过程中，保存每个检查点时刻及前一时刻的p和q全部波场。参见图2中深灰色、浅灰色和白色区域为p和q全部波场将其存储至计算机内。比如，当震源波场正传至检查点N时，记录tNc-1和tNc两个时刻的p和q全部波场。如图3(a)所示。所述每个检查点时刻及前一时刻的p和q全部波场作为存储若干部分波场进行存储。
优选地，所述步骤S4，具体为：
先在所述若干部分波场中确定N个检查点；其中若干部分波场和检查点的个数是相等的。
在对所述检查点N-1至检查点N时段内波场进行重建时，读取所述检查点N-1时刻及前一时刻的p和q全部波场信息作为初值进行c次正向延拓，每次延拓记录内边界层信息，保存至图形处理器的显存p和q缓冲区，所述显存p和q缓冲区为两个检查点之间p和q内边界大小；此处，图形处理器的显存存储空间只是两个检查点之间p和q内边界的大小，并不是所有时间p和q内边界的大小。也就是说对每两个检查点间实现波场重建的过程中，只需更新图形处理器显存中的p和q波场的内边界信息即可。
将tNc和tNc-1两时刻的p和q全部波场作为反向延拓的初值，每延拓一次，将所述保存至显存上p和q缓冲区对应时刻内边界信息载入重建波场中，反向延拓c-1次得到重建t(N-1)c时刻波场；
重复上述步骤重建tNc至t0波场，得到震源波场。
其中，c的选取基于存储空间进行选取，如果c的数值大，检查点数就少，数据交换的次数就少，CPU占用内存就小，但是GPU显存占用就大。
受有限差分数值解法的限制，在计算内边界层(如图2浅灰色区域)波场信息时，需要利用完全匹配吸收边界的信息，由于边界内信息是吸收衰减后的结果，因此逆向延拓得到的tNc-2时刻的p和q内边界层波场(如图2浅灰色区域)与带重建波场存在误差，因为在逆时偏移成像过程中只需保证Nx×Nz区域波场精度即可。
优选地，在进行所述正向延拓和逆向延拓时，通过图形处理器进行计算。
由于求解qP波方程计算量大，本实施例中采用GPU加快计算速度。
另外，在本实施例中，举例说明对qP波波动方程进行数值模拟，及经过本发明方法进行模拟的效果说明。
设定一个模拟模型，假定模型在水平和垂直方向的网格数均为300，水平方向和垂直方向的网格间距均为10m，沿对称轴方向的纵波速度Vpz＝3000m/s，完全匹配吸收层为50，震源位于模型的中央，主频为30Hz的雷克子波做为激发震源，时间采样间隔为0.5ms。图4(a)为300ms的波场图像，其中两个各向异性参数为零，则公式(1)退化为声波方程；图4(b)中横波的速度为0，在波场快照中可以观察到菱形尖角状的伪横波，在各向异性参数变化剧烈的情况下，该伪横波将导致严重的频散；图4(c)为在公式(1)中引入了横波速度，菱形尖角状消失；图4(d)为600ms的波场图像，证明所加的吸收匹配层正确有效，消除自由边界的反射现象。
其中，图4(a)为vsz＝0，ε＝0，δ＝0，t＝300ms时的波场图像，图4(b)为vsz＝0，ε＝0.3，δ＝0.1，t＝300ms时的波场图像，图4(c)为vsz＝1549m/s，ε＝0.3，δ＝0.1，t＝300ms时的波场图像；图4(d)为vsz＝1549m/s，ε＝0.3，δ＝0.1，t＝600ms时的波场图像。
为了验证检查点存储策略的正确有效性，设计VTI介质模型，速度如图5所示，将震源位于模型中间位置，现得到1150ms的波场图像，如图6(a)所示，图6(b)为利用本发明中的检查点波场存储策略对1150ms的波场实现重建的结果，将图6(a)和图6(b)做差得到图6(c)，该图中的振幅值在10-5，可以认为是系统误差，重建精度较高，检查点存储策略更加适用于大数据模型在GPU端的计算。
通过本发明的方法对复杂VTI介质模型进行逆时偏移成像计算，图7(a)为介质的速度，图7(b)和图7(c)分别为各向异性参数δ和ε，图7(d)为逆时偏移成像剖面，图中可见：模型中构造体空间归位准确，断层边缘清晰，实现了对VTI介质的高精度成像。
第二方面，本发明提供一种针对VTI介质逆时偏移设备，包括：至少一个处理器、至少一个存储器以及存储在所述存储器中的计算机程序指令，当所述计算机程序指令被所述处理器执行时实现如第一方面所述的方法。
结合图8描述的本发明实施例的针对VTI介质逆时偏移方法可以由针对VTI介质逆时偏移设备来实现。图8示出了本发明实施例提供的针对VTI介质逆时偏移设备的硬件结构示意图。
针对VTI介质逆时偏移设备可以包括处理器401以及存储有计算机程序指令的存储器402。
具体地，上述处理器401可以包括中央处理器(CPU)，或者特定集成电路(Application Specific Integrated Circuit，ASIC)，或者可以被配置成实施本发明实施例的一个或多个集成电路。
存储器402可以包括用于数据或指令的大容量存储器。举例来说而非限制，存储器402可包括硬盘驱动器(Hard Disk Drive，HDD)、软盘驱动器、闪存、光盘、磁光盘、磁带或通用串行总线(Universal Serial Bus，USB)驱动器或者两个或更多个以上这些的组合。在合适的情况下，存储器402可包括可移除或不可移除(或固定)的介质。在合适的情况下，存储器402可在数据处理装置的内部或外部。在特定实施例中，存储器402是非易失性固态存储器。在特定实施例中，存储器402包括只读存储器(ROM)。在合适的情况下，该ROM可以是掩模编程的ROM、可编程ROM(PROM)、可擦除PROM(EPROM)、电可擦除PROM(EEPROM)、电可改写ROM(EAROM)或闪存或者两个或更多个以上这些的组合。
处理器401通过读取并执行存储器402中存储的计算机程序指令，以实现上述实施例中的任意一种针对VTI介质逆时偏移方法。
在一个示例中，针对VTI介质逆时偏移设备还可包括通信接口403和总线410。其中，如图8所示，处理器401、存储器402、通信接口403通过总线410连接并完成相互间的通信。
通信接口403，主要用于实现本发明实施例中各模块、装置、单元和/或设备之间的通信。
总线410包括硬件、软件或两者，将针对VTI介质逆时偏移设备的部件彼此耦接在一起。举例来说而非限制，总线可包括加速图形端口(AGP)或其他图形总线、增强工业标准架构(EISA)总线、前端总线(FSB)、超传输(HT)互连、工业标准架构(ISA)总线、无限带宽互连、低引脚数(LPC)总线、存储器总线、微信道架构(MCA)总线、外围组件互连(PCI)总线、PCI-Express(PCI-X)总线、串行高级技术附件(SATA)总线、视频电子标准协会局部(VLB)总线或其他合适的总线或者两个或更多个以上这些的组合。在合适的情况下，总线410可包括一个或多个总线。尽管本发明实施例描述和示出了特定的总线，但本发明考虑任何合适的总线或互连。
第三方面，结合上述实施例中的针对VTI介质逆时偏移方法，本发明实施例可提供一种计算机可读存储介质来实现。该计算机可读存储介质上存储有计算机程序指令；该计算机程序指令被处理器执行时实现上述实施例中的任意一种针对VTI介质逆时偏移方法。
需要明确的是，本发明并不局限于上文所描述并在图中示出的特定配置和处理。为了简明起见，这里省略了对已知方法的详细描述。在上述实施例中，描述和示出了若干具体的步骤作为示例。但是，本发明的方法过程并不限于所描述和示出的具体步骤，本领域的技术人员可以在领会本发明的精神后，作出各种改变、修改和添加，或者改变步骤之间的顺序。
以上所述的结构框图中所示的功能块可以实现为硬件、软件、固件或者它们的组合。当以硬件方式实现时，其可以例如是电子电路、专用集成电路(ASIC)、适当的固件、插件、功能卡等等。当以软件方式实现时，本发明的元素是被用于执行所需任务的程序或者代码段。程序或者代码段可以存储在机器可读介质中，或者通过载波中携带的数据信号在传输介质或者通信链路上传送。“机器可读介质”可以包括能够存储或传输信息的任何介质。机器可读介质的例子包括电子电路、半导体存储器设备、ROM、闪存、可擦除ROM(EROM)、软盘、CD-ROM、光盘、硬盘、光纤介质、射频(RF)链路，等等。代码段可以经由诸如因特网、内联网等的计算机网络被下载。
还需要说明的是，本发明中提及的示例性实施例，基于一系列的步骤或者装置描述一些方法或系统。但是，本发明不局限于上述步骤的顺序，也就是说，可以按照实施例中提及的顺序执行步骤，也可以不同于实施例中的顺序，或者若干步骤同时执行。
与现有技术相比，本发明的优点：
1、常规声波方程只能对地震波在各向同性介质进行精确描述，本发明利用准纵波方程实现地震波在VTI介质中的数值模拟，与利用弹性波方程相比准纵波方程在大幅减小计算量的情况下，能够获得相同运动学特征的纵波分量。
2、在准纵波方程中加入稳定性因子，该因子能够消除数值模拟中波场出现的频散现象，确保地震波在剧烈变化的VTI介质中能够稳定地延拓，为获取高精度逆时偏移成像剖面提供前提保障。
3、对于波场的有限差分数值解法，从时间和空间两个角度建立一种检查点存储策略，从而能够极大地节约波场存储空间，同时该策略能降低传统波场重建带来的误差积累，降低逆时偏移算法对计算机硬件的依赖，为方法的工业化带来更广阔的空间。
4、将图形处理器(GPU)设备应用到波场延拓计算中，利用中央处理器设备逻辑判断能力强和GPU设备的计算效率高的特点，形成CPU/GPU协同加速技术，从而大幅度减小计算周期。
5、基于qP波方程建立针对VTI介质逆时偏移成像的处理流程，使成像剖面的反射同相轴准确归位，提升复杂介质的成像精度。
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。
Claims (6)
1.一种针对VTI介质逆时偏移方法，其特征在于，包括：
步骤S1，获取炮集数据；
步骤S2，在所述炮集数据中提取震源坐标，并生成震源子波；
步骤S3，根据正演方案，将震源子波加载至波场中，得到震源波场，通过qP波波动方程对所述震源波场沿时间正向延拓，共进行Nc次延拓，Nc为t0至tNc时间段内设定的延拓次数，并在延拓过程中选择若干部分波场进行存储；
在所述qP波波动方程中引入稳定性因子σ；
利用稳定性因子公式实现对横波沿着对称轴方向的传播速度Vsz的计算，其中，通过数值模拟，所述稳定性因子σ为0.75，ε和δ为Thomsen参数，Vpz为纵波沿着对称轴方向的传播速度；
所述步骤S3中，根据正演方案，通过qP波波动方程对所述震源子波沿时间正向延拓，具体为：
根据所述纵波和横波沿着对称轴方向的传播速度以及Thomsen参数作为基于qP波波动方程的VTI介质逆时偏移成像方法的模型参数；
将地震子波作为激发震源生成地震波场，利用有限差分数值解法，对所述震源波场沿时间正向Nc次延拓，Nc为t0至tNc时间段内设定的延拓次数；
其中，qP波波动方程如以下公式所示，
式中，Vpz和Vsz分别为纵波和横波沿着对称轴方向的传播速度，ε和δ为Thomsen参数，Vpx为纵波沿着垂直于对称轴方向的传播速度，Vpn为动校正速度，p为声压，q为辅助变量；
步骤S4，根据所述正演方案，通过所述qP波波动方程对所述若干部分波场沿时间逆向延拓，共进行Nc次延拓，实现对所述震源波场的重建；
假设所述震源波场的正向延拓的总时间为tNc，在所述震源波场中确定N个检查点；
在所述震源波场正向延拓的过程中，保存每个检查点时刻及前一时刻的p和q全部波场；
将所述每个检查点时刻及前一时刻的p和q全部波场作为若干部分波场进行存储，其中，p为声压，q为辅助变量；
所述步骤S4，具体为：
在对所述检查点N-1至检查点N时段内波场进行重建时，读取所述检查点N-1时刻及前一时刻的p和q全部波场信息作为初值进行c次正向延拓，每次延拓记录内边界层信息，保存至图形处理器的显存p和q缓冲区，所述显存p和q缓冲区为两个检查点之间p和q内边界大小；
将tNc和tNc-1两时刻的p和q全部波场作为反向延拓的初值，每延拓一次将所述保存至显存上p和q缓冲区对应时刻内边界信息载入重建波场中，反向延拓c-1次得到重建t(N-1)c时刻波场；
重复上述步骤重建tNc至t0波场，得到震源波场；
步骤S5，在所述炮集数据中提取地震数据，并对所述地震数据进行插值，得到插值后的地震数据；
步骤S6，根据所述正演方案，通过所述qP波波动方程对所述插值后的地震数据沿时间逆向延拓，共进行Nc次延拓，得到检波点波场；
步骤S7，将所述震源波场和所述检波点波场进行互相关成像，得到一炮偏移剖面，重复上述步骤，得到多炮偏移剖面；
步骤S8，将所述多炮偏移剖面进行叠加，得到最终成像剖面。
2.根据权利要求1所述的方法，其特征在于，
所述步骤S2之后，还包括：
对所述震源子波进行插值，得到插值后的震源子波。
3.根据权利要求1所述的方法，其特征在于，
基于qP波波动方程进行数值模拟计算具体为：
利用高阶有限差分数值方法求解qP波波动方程，采用PML完全匹配层边界条件吸收边界反射，得到模拟地震波场。
4.根据权利要求1～3任一项所述的方法，其特征在于，
在进行所述正向延拓和逆向延拓时，通过图形处理器进行计算。
5.一种针对VTI介质逆时偏移设备，其特征在于，包括：至少一个处理器、至少一个存储器以及存储在所述存储器中的计算机程序指令，当所述计算机程序指令被所述处理器执行时实现如权利要求1～3中任一项所述的方法。
6.一种计算机可读存储介质，其上存储有计算机程序指令，其特征在于，当所述计算机程序指令被处理器执行时实现如权利要求1-3中任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810039657.7A CN108181653B (zh) | 2018-01-16 | 2018-01-16 | 针对vti介质逆时偏移方法、设备及介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810039657.7A CN108181653B (zh) | 2018-01-16 | 2018-01-16 | 针对vti介质逆时偏移方法、设备及介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108181653A CN108181653A (zh) | 2018-06-19 |
| CN108181653B true CN108181653B (zh) | 2019-11-19 |
Family
ID=62550601
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810039657.7A Expired - Fee Related CN108181653B (zh) | 2018-01-16 | 2018-01-16 | 针对vti介质逆时偏移方法、设备及介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108181653B (zh) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109490955B (zh) * | 2018-11-14 | 2021-07-20 | 深圳市勘察研究院有限公司 | 一种基于规则网格的声波波动方程正演模拟方法及装置 |
| CN109901221B (zh) * | 2019-03-05 | 2021-03-19 | 中国石油大学(华东) | 一种基于动校正速度参数的地震资料各向异性建模方法 |
| CN110618459B (zh) * | 2019-11-01 | 2020-11-10 | 中国科学院地质与地球物理研究所 | 地震数据处理方法和装置 |
| CN111487677A (zh) * | 2020-03-31 | 2020-08-04 | 深圳市勘察研究院有限公司 | 一种声波波动方程叠前逆时偏移成像方法及装置 |
| CN112213783B (zh) * | 2020-09-16 | 2021-11-12 | 河海大学 | 一种利用海底多分量地震记录进行各向异性地震成像方法 |
| CN112987088B (zh) * | 2021-02-22 | 2023-04-18 | 成都理工大学 | 一种渗流介质地震横波数值模拟和成像方法 |
| CN117233838B (zh) * | 2023-09-20 | 2024-04-05 | 长江大学 | 一种二维vti介质中的弹性准纵横波场分离和逆时偏移成像方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102590859A (zh) * | 2011-12-31 | 2012-07-18 | 中国石油集团西北地质研究所 | Vti介质准p波方程各向异性逆时偏移方法 |
| CN105242313A (zh) * | 2015-09-06 | 2016-01-13 | 中国科学院地质与地球物理研究所 | 一种弹性波逆时偏移极性反转的校正方法及系统 |
| CN106597535A (zh) * | 2016-12-01 | 2017-04-26 | 桂林理工大学 | 一种提高弹性波逆时偏移计算率和空间分辨率的方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101903475B1 (ko) * | 2016-02-04 | 2018-10-02 | 한국해양대학교 산학협력단 | 수직횡적등방성 매질의 암석물성정보 분석 방법 및 장치 |
-
2018
- 2018-01-16 CN CN201810039657.7A patent/CN108181653B/zh not_active Expired - Fee Related
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102590859A (zh) * | 2011-12-31 | 2012-07-18 | 中国石油集团西北地质研究所 | Vti介质准p波方程各向异性逆时偏移方法 |
| CN105242313A (zh) * | 2015-09-06 | 2016-01-13 | 中国科学院地质与地球物理研究所 | 一种弹性波逆时偏移极性反转的校正方法及系统 |
| CN106597535A (zh) * | 2016-12-01 | 2017-04-26 | 桂林理工大学 | 一种提高弹性波逆时偏移计算率和空间分辨率的方法 |
Non-Patent Citations (2)
| Title |
|---|
| VTI介质高精度叠前逆时偏移方法研究;方修改 等;《地球物理学报》;20151231;第30卷(第4期);1677-1684 * |
| 基于一步法波场延拓的正演模拟和逆时偏移成像;柯璇 等;《地球物理学报》;20171130;第60卷(第11期);4468-4479 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108181653A (zh) | 2018-06-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108181653B (zh) | 针对vti介质逆时偏移方法、设备及介质 | |
| CN104297785B (zh) | 岩相约束储层物性参数反演方法及装置 | |
| CN105137486B (zh) | 各向异性介质中弹性波逆时偏移成像方法及其装置 | |
| CN103293552B (zh) | 一种叠前地震资料的反演方法及系统 | |
| US20170293041A1 (en) | Device, system and method for geological-time refinement | |
| CN105425289B (zh) | 确定低频波阻抗的方法和装置 | |
| CN110082823B (zh) | 一种地震数据插值方法及装置 | |
| CN113064203B (zh) | 共轭梯度归一化lsrtm方法、系统、存储介质及应用 | |
| CN103149586A (zh) | 一种倾斜层状粘弹性介质中波场正演模拟方法 | |
| CN111596366A (zh) | 一种基于地震信号优化处理的波阻抗反演方法 | |
| CN110058303A (zh) | 声波各向异性逆时偏移混合方法 | |
| CN107765308A (zh) | 基于褶积思想与精确震源的重构低频数据频域全波形反演方法 | |
| CN109946742A (zh) | 一种TTI介质中纯qP波地震数据模拟方法 | |
| CN112231974A (zh) | 基于深度学习的tbm破岩震源地震波场特征恢复方法及系统 | |
| CN106353798A (zh) | 多分量联合高斯束叠前逆时偏移成像方法 | |
| CN110515122B (zh) | 正演网格搜索定位及微地震信号识别方法及装置 | |
| CN111582114A (zh) | 一种地震断层识别方法、装置、设备和存储介质 | |
| CN106574980A (zh) | 用于地下地质体的岩石性质估计的系统和方法 | |
| CN105353409A (zh) | 一种用于抑制全波形反演震源编码串扰噪音的方法和系统 | |
| CN106646593B (zh) | 一种跨节点并行的三维起伏地表声波正演模拟方法 | |
| US20160034612A1 (en) | Re-ordered Interpolation and Convolution for Faster Staggered-Grid Processing | |
| Zhou et al. | An iterative factored topography-dependent eikonal solver for anisotropic media | |
| CN115327626A (zh) | 一种用于三维ati介质的弹性波场矢量分解方法及系统 | |
| CN115600373A (zh) | 黏滞各向异性介质qP波模拟方法、系统、设备及应用 | |
| CN114895351A (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 | ||
| 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: 20191119 Termination date: 20220116 |