CN116413801B - 一种各向异性介质弹性波高精度成像方法和系统 - Google Patents
一种各向异性介质弹性波高精度成像方法和系统 Download PDFInfo
- Publication number
- CN116413801B CN116413801B CN202310105843.7A CN202310105843A CN116413801B CN 116413801 B CN116413801 B CN 116413801B CN 202310105843 A CN202310105843 A CN 202310105843A CN 116413801 B CN116413801 B CN 116413801B
- Authority
- CN
- China
- Prior art keywords
- wave
- wave field
- anisotropic
- vector
- elastic wave
- 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
- 238000003384 imaging method Methods 0.000 title claims abstract description 58
- 230000002441 reversible effect Effects 0.000 claims abstract description 61
- 238000013508 migration Methods 0.000 claims abstract description 55
- 230000005012 migration Effects 0.000 claims abstract description 48
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 38
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 16
- 230000002194 synthesizing effect Effects 0.000 claims abstract description 13
- 230000006870 function Effects 0.000 claims description 87
- 238000000034 method Methods 0.000 claims description 36
- 238000004364 calculation method Methods 0.000 claims description 21
- 238000004088 simulation Methods 0.000 claims description 12
- 239000011159 matrix material Substances 0.000 claims description 6
- 238000002939 conjugate gradient method Methods 0.000 claims description 3
- 238000001514 detection method Methods 0.000 claims description 3
- 238000001615 p wave Methods 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 38
- 238000006073 displacement reaction Methods 0.000 description 13
- 239000002245 particle Substances 0.000 description 11
- 238000011161 development Methods 0.000 description 9
- 238000004590 computer program Methods 0.000 description 7
- 238000005070 sampling Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 4
- 238000004519 manufacturing process Methods 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 239000004215 Carbon black (E152) Substances 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 229930195733 hydrocarbon Natural products 0.000 description 3
- 150000002430 hydrocarbons Chemical class 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000005286 illumination Methods 0.000 description 2
- 230000010287 polarization Effects 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 229920002595 Dielectric elastomer Polymers 0.000 description 1
- 206010027336 Menstruation delayed Diseases 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 150000003839 salts Chemical class 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 230000001568 sexual effect Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 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/30—Analysis
- G01V1/301—Analysis for determining seismic cross-sections or geostructures
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
本发明涉及一种各向异性介质弹性波高精度成像方法和系统：求解各向异性弹性波动方程，得到震源端波场；对震源端波场进行矢量分解，获得P波矢量波场和S波矢量波场；进行各向异性弹性波逆时反偏移，合成多分量反射地震数据；计算多分量反射地震数据与多分量实际地震数据的目标函数和残差，若其满足收敛条件输出高分辨率图像，否则进入下一步；将残差进行各向异性弹性波逆时偏移，获得目标函数的梯度；计算步长和共轭梯度方向，结合梯度得到共轭梯度算法，更新多分量弹性波图像模型，重新合成多分量反射地震数据，重复上述步骤直至目标函数满足收敛条件。其获得的各向异性介质弹性波PP、PS、SP、SS图像剖面，为复杂油气藏提供更完备的成像信息。
Description
技术领域
本发明涉及一种各向异性介质弹性波高精度成像方法和系统，属于地质勘探技术领域，特别涉及地震偏移成像技术领域。
背景技术
地震偏移成像利用布置在地表或井中的接收器观测的地震数据来建立地下介质反射系数图像。地震偏移成像对油气藏的勘探开发具有重要意义。在油气勘探阶段，地震偏移成像是探测油气藏位置主要手段之一。在油气开发阶段，地震偏移图像可用来反演地球介质弹性参数信息，指导后续叠前储层反演和油气藏地质建模。在油气生产后期，基于时移地震采集技术的地震成像可用于油气藏的动态监测来提高油气采收率等。经过长期的勘探开发，我国易于开发的高品质油气藏逐渐减少，复杂地质环境的油气藏资源的勘探开发对提高我国油气资源自给能力和保障国家能源安全至关重要。
传统的地震成像技术基于声波假设，地球介质被近似为流体，仅由一个参数即纵波速度来描述。然而，地球介质是弹性介质，其中纵横波的传播是相互耦合的矢量波场。地震横波是描述油气储层特征的有用工具。例如，比照PP和PS反射系数图像能帮助识别油水界面，转换横波射线路径能提供气烟囱下的储层照明，横波连同纵波能提供完整的数据信息来反演弹性参数。此外，相同观测采集系统下，转换横波比纵波有更宽的源端照明成像孔径，横波图像比纵波图像具有更高的空间分辨率等。近年来，弹性波逆时偏移成像得到了一些发展，但主要基于各向同性的弹性介质假设，未考虑介质的各向异性特性。
复杂油气藏具有非均质性强、构造复杂和弹性各向异性特性，对地震偏移成像带来了巨大的挑战。地下介质的各向异性会导致地震波的走时和振幅的变化。如果不考虑介质的各向异性，地震偏移图像存在反射点深度偏差和图像振幅响应的不准确性，影响对后续的地震解释，从而造成油气藏的定位偏差。纵波(P波)和横波(S波)波场分离是各向异性弹性波逆时偏移的重要步骤。然而，当弹性介质存在各向异性特性时，纵横波的偏振方向既不平行也不垂直于它们的传播方向，偏振方向、群速度方向和相速度方向三者可能皆互不相同，各向同性亥姆霍兹分解方法不再能准确地解耦P/S波模式。
近年来，基于弹性波场矢量分解的各向异性弹性波逆时偏移成像技术成为科研机构和石油天然气公司的研发前沿。然而，各向异性弹性波逆时偏移只是正演算子的伴随算子，并不是其逆算子。故现有技术中大部分都是基于波动方程的线性化正演算子的伴随过程。基于正演算子的伴随算子的各向异性弹性波逆时偏移存在如下问题和挑战：(1)偏移图像存在空间分辨率相对较低、低频噪声和波模式串扰等问题；(2)当数据存在噪声，或因震源检波器采样差导致数据缺失情况时，偏移成像的效果较差；(3)偏移图像振幅不均衡，不能准确地反映反射系数，影响后续角道集AVA(Amplitude Versus Angle)分析。如今，复杂油气藏的勘探开发面临巨大的挑战。油气工业界现有的地震成像技术难以满足当前复杂油气藏的高精度勘探开发的需求。
发明内容
针对上述问题，本发明的目的是提供一种各向异性介质弹性波高精度成像方法和系统，其基于弹性波场矢量分解的各向异性介质弹性波最小二乘逆时偏移成像，利用多分量地震勘探数据得到各向异性介质高精度多分量弹性波PP、PS、SP、SS图像剖面，为复杂油气藏提供更完备的成像信息，能更好地进行油气储层描述和地震解释。
为实现上述目的，本发明提出了以下技术方案：一种各向异性介质弹性波高精度成像方法，包括以下步骤：通过交错网格有限差分，根据各向异性参数模型，求解各向异性弹性波动方程，得到震源端波场；对所述震源端波场进行矢量分解，获得P波矢量波场和S波矢量波场；利用震源端的所述P波矢量波场、S波矢量波场和输入的多分量弹性波图像初始模型，进行各向异性弹性波逆时反偏移，合成多分量反射地震数据；计算所述多分量反射地震数据与多分量实际地震数据的目标函数，若所述目标函数满足收敛条件输出高分辨率图像，若所述目标函数不满足收敛条件，则进入下一步；计算所述多分量反射地震数据与多分量实际地震数据的残差，将所述残差进行各向异性弹性波逆时偏移，获得目标函数的梯度；计算步长和共轭梯度方向，结合所述梯度得到共轭梯度算法，采用所述共轭梯度算法更新多分量弹性波图像模型，重新合成多分量反射地震数据，计算目标函数，并判断目标函数是否满足收敛条件，直至目标函数满足收敛条件。
进一步，所述目标函数为：
其中，m(i)是当前第(i)次迭代的多分量弹性波图像；L是各向异性弹性波逆时反偏移算子，dobs是多炮多分量地震数据矢量；如果目标函数小于设定的阈值，则输出高分辨率图像，如果目标函数大于等于设定的阈值，则进入下一步。
进一步，所述震源端波场通过计算以下异性弹性波动方程获得：
其中，u＝(ux,uz)T是震源端波场，f＝(fx,fz)T是矢量震源，ρ是密度，Cij表示刚度矩阵的元素。
进一步，对所述震源端波场进行矢量分解的方法为：通过快速傅里叶变换把弹性波场时间切片变换到波数-空间混合域，并求解各向异性矢量泊松方程，得到辅助矢量波场w，其计算公式为：
P波矢量波场uP(x，z)的计算公式为：
S波矢量波场uS(x，z)的计算公式为：
其中，是各向异性空间梯度算子。
进一步，由各向异性弹性波逆时反偏移获得所述多分量反射地震数据的方法为：将震源端波场进行矢量分解后获得的P波矢量波场和S波矢量波场与对应的波模式反射系数图像相乘，得到波场-反射系数图像乘积；将所述波场-反射系数图像乘积分解成四种波模式，得到反射点处的PP波场，PS波场，SP波场和SS波场；将所述四种波模式结合作为二次震源，对所述二次震源进行基于交错网格有限差分的正向传播数值模拟，在检波点处抽取由反偏移算子合成的反射波，合成多分量反射地震数据。
进一步，所述各向异性弹性波逆时偏移的具体方法为：将所述残差作为伴随震源，通过波场反向传播模拟，得到检波器端波场；通过弹性波场矢量分解，将所述检波器端波场分解成P波矢量波场和S波矢量波场；由所述震源端波场和检波器端波场，根据矢量成像条件把对应的波模互相关生成PP、PS、SP、SS偏移图像，并生成目标函数的梯度。
进一步，所述检波器端波通过计算以下异性弹性波动方程获得：
其中，φ＝(φx，φz)T是检波器端波场，r＝(rx，rz)T是伴随震源，ρ是密度，Cij表示刚度矩阵的元素。
进一步，所述目标函数的梯度的公式为：
其中，g＝(gPP，gPS，gSP，gSS)T是目标函数的梯度，uP和φP分别是解耦后的震源端、检波器端P波矢量波场，uS和φS分别是解耦后的震源端、检波器端S波矢量波场；Tmax是地震数据的最大接收时间。
进一步，所述步长α的计算公式为：
共轭梯度方向c(i+1)的计算公式为：
c(i+1)＝Pg(i+1)+βc(i)
其中，L是各向异性弹性波逆时反偏移算子，P是预条件算子，β是系数。
本发明还公开了一种各向异性介质弹性波高精度成像系统，包括：波场正向传播模拟模块，用于通过交错网格有限差分，根据各向异性参数模型，求解各向异性弹性波动方程，得到震源端波场；弹性波场矢量分解模块，用于对所述震源端波场进行矢量分解，获得P波矢量波场和S波矢量波场；各向异性弹性波逆时反偏移模块，用于通过震源端的所述P波矢量波场、S波矢量波场和输入的多分量弹性波图像初始模型，进行各向异性弹性波逆时反偏移，合成多分量反射地震数据；目标函数判断模块，用于计算所述多分量反射地震数据与多分量实际地震数据的目标函数，若所述目标函数满足收敛条件输出高分辨率图像，若所述目标函数不满足收敛条件，则进入下一模块；目标函数梯度计算模块，用于计算所述多分量反射地震数据与多分量实际地震数据的残差，将所述残差进行各向异性弹性波逆时偏移，获得目标函数的梯度；共轭梯度法偏移图像更新模块，用于计算步长和共轭梯度方向，结合所述梯度得到共轭梯度算法，采用所述共轭梯度算法更新多分量弹性波图像模型，重新合成多分量反射地震数据，计算目标函数，并判断目标函数是否满足收敛条件，直至目标函数满足收敛条件。
本发明由于采取以上技术方案，其具有以下优点：
1、本发明能够利用多分量地震勘探数据得到各向异性介质高精度多分量弹性波PP、PS、SP、SS图像剖面，为复杂油气藏提供更完备的成像信息，能更好地进行油气储层描述和地震解释；
2、本发明中方案把反演理论引入各向异性介质弹性波成像中，基于弹性波场矢量分解的各向异性弹性波最小二乘逆时偏移方法；相比传统的各向异性弹性波逆时偏移，本发明提高了图像的空间分辨率，减少了低频逆时偏移噪声，压制了波模式串扰噪声，显著提高了成像质量；
3、本发明的基于弹性波场矢量分解的各向异性弹性波最小二乘逆时偏移方法，基于最小二乘反演理论，在地震数据存在噪声或采样缺失时，能够嵌入数据权重函数，对质量差的数据赋予低的权重，减少对图像的影响。
4、本发明中方案为复杂油气藏精细勘探和开发提供支撑，助力对复杂油气藏的开发和利用，实现增储上产。
附图说明
图1是本发明一实施例中各向异性介质弹性波高精度成像方法的流程图；
图2是本发明一实施例中各向异性介质弹性波逆时反偏移模块的流程示意图；
图3是本发明一实施例中目标函数梯度计算模块的流程示意图；
图4是Hess各向异性弹性模型的示意图，图4(a)是纵波速度模型示意图，图4(b)是横波速度模型示意图，图4(c)是Thomsen各向异性参数ε模型示意图，图4(d)为Thomsen各向异性参数δ模型示意图，图4(e)为密度模型示意图；
图5是本发明一实施例中多分量实际地震数据，图5(a)是水平分量位移波场数据，图5(b)是垂直分量位移波场数据；
图6是本发明一实施例中基于各向异性亥姆霍兹分解方法对震源端弹性波场进行波场矢量分解图；图6(a)是水平分量质点位移波场图，图6(b)是垂直分量质点位移波场图，图6(c)是解耦后的水平分量P波场图，图6(d)是解耦后的垂直分量P波场图，图6(e)是解耦后的水平分量S波场图，图6(f)为解耦后的垂直分量S波场图；
图7是本发明一实施例中基于各向异性亥姆霍兹分解方法对检波器端弹性波场进行波场矢量分解图；图7(a)是水平分量质点位移波场图，图7(b)是垂直分量质点位移波场图，图7(c)是解耦后的水平分量P波场图，图7(d)是解耦后的垂直分量P波场图，图7(e)是解耦后的水平分量S波场图，图7(f)为解耦后的垂直分量S波场图；
图8是现有技术中VTI(Vertical transverselyisotropic)各向异性弹性波逆时偏移成像图；图8(a)是PP图，图8(b)是PS图，图8(c)是SP图，图8(d)是SS图；
图9是本发明一实施例中VTI各向异性弹性波最小二乘逆时偏移成像图；图9(a)是PP图，图9(b)是PS图，图9(c)是SP图，图9(d)是SS图；
图10是本发明一实施例中VTI各向异性弹性波最小二乘逆时偏移的数据误差收敛曲线示意图。
具体实施方式
为了使本领域技术人员更好的理解本发明的技术方案，通过具体实施例对本发明进行详细的描绘。然而应当理解，具体实施方式的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。在本发明的描述中，需要理解的是，所用到的术语仅仅是用于描述的目的，而不能理解为指示或暗示相对重要性。
为了解决现有技术中存在的偏移图像存在空间分辨率相对较低、低频噪声和波模式串扰等问题；当数据存在噪声，或因震源检波器采样差导致数据缺失情况时，偏移成像的效果较差；偏移图像振幅不均衡，不能准确地反映反射系数，影响后续角道集AVA分析。本发明提出了各向异性介质弹性波高精度成像方法和系统；其基于弹性波场矢量分解的逆时偏移和反偏移算子，利用多分量地震勘探数据得到各向异性介质高精度多分量弹性波PP、PS、SP、SS图像剖面，为复杂油气藏提供更完备的成像信息，能更好地进行油气储层描述和地震解释。其把反演理论引入各向异性介质弹性波成像中，发明了基于弹性波场矢量分解的各向异性弹性波最小二乘逆时偏移方法；相比传统的各向异性弹性波逆时偏移，其提高了图像的空间分辨率，减少了低频逆时偏移噪声，压制了波模式串扰噪声，大大提高了成像质量；基于最小二乘反演理论，在地震数据存在噪声或采样缺失时，能够嵌入数据权重函数，对质量差的数据赋予低的权重，减少对图像的影响。下面结合附图，通过实施例对本发明方案进行详细阐述。
实施例一
本实施例公开了一种各向异性介质弹性波高精度成像方法，如图1所示，包括以下步骤：
S1输入地层密度、垂向纵波横波速度、Thomsen各向异性参数的偏移背景模型和多分量弹性波图像(PP、PS、SP、SS)的初始值m(0)，通过交错网格有限差分求解各向异性弹性波动方程，得到震源端波场，如图2所示；
震源端波场通过计算以下各向异性弹性波动方程获得：
其中，u＝(ux，uz)T是震源端波场，f＝(fx，fz)T是矢量震源，ρ是密度，Cij表示刚度矩阵的元素，刚度矩阵模型可由输入的垂向纵横波速度和Thomsen各向异性参数转换得到。
S2基于各向异性亥姆霍兹分解方法，对震源端波场进行矢量分解，获得P波矢量波场和S波矢量波场；
对震源端波场进行矢量分解的方法为：通过基于快速傅里叶变换的波数-空间混合域方法，求解各向异性矢量泊松方程，得到辅助矢量波场w，其计算公式为：
P波矢量波场uP(x，z)的计算公式为：
S波矢量波场uS(x，z)的计算公式为：
其中，是各向异性空间梯度算子，由克里斯托弗尔方程的特征向量推导得到。
其中，和/>分别是对x和z方向的偏导，r是系数，∈和δ是Thomsen各向异性参数，Vp0和Vs0分别是垂向纵波和横波速度，/>是各向异性散度；/>是各向异性旋度算子。
获得多分量反射地震数据的方法为：
S2.1获得的P波矢量波场和S波矢量波场与对应的波模式反射系数图像相乘(mPP，mPS，mSP，mSS)，得到波场-反射系数图像乘积；
S2.2将波场-反射系数图像乘积分解成四种波模式，得到反射点处的PP波场，PS波场，SP波场和SS波场；
S2.3将四种波模式结合作为二次震源，对二次震源进行基于交错网格有限差分的正向传播数值模拟，即以二次震源作为波动方程的震源函数，求解如下波场正向模拟方程；在检波点处抽取由反偏移算子合成的反射波，合成多分量反射地震数据，由公式表示为dsyn＝Lm(i)，其中，L代表各向异性弹性波逆时反偏移算子，m(i)表示第(i)次迭代的多分量弹性波图像。
其中，δu＝(δux，δuz)T是散射的质点位移矢量波场，s＝(sx，sz)T表示二次震源。
S3从电脑磁盘读入SEGY格式多分量实际地震数据炮集记录dobs，计算多分量反射地震数据与多分量实际地震数据dobs的残差和目标函数，若目标函数满足收敛条件输出高分辨率图像，若目标函数不满足收敛条件，则进入下一步；
可选的，本实施例中的目标函数为：
其中，m(i)是当前第(i)次迭代的多分量弹性波图像；L是各向异性弹性波逆时反偏移算子，dobs是多炮多分量地震数据矢量；如果目标函数小于设定的阈值，则输出高分辨率图像，如果目标函数大于等于设定的阈值，则进入下一步。
S4计算多分量反射地震数据与多分量实际地震数据的残差r(i)＝dobs-Lm(i)，根据数据残差获得检波器端波场，根据震源端波场和检波器端波场获得目标函数的梯度；
由各向异性弹性波逆时偏移得到目标函数的梯度的具体方法，如图3所示，为：
S4.1将残差作为伴随震源，通过波场反向传播模拟，即通过交错网格有限差分方法求解如下伴随波动方程得到检波器端波场φ；
其中，φ＝(φx，φz)T是检波器端波场，r＝(rx，rz)T是伴随震源，ρ是密度，Cij表示刚度矩阵的元素。
S4.2通过弹性波场矢量分解，将检波器端波场φ分解成P波矢量波场φP和S波矢量波场φS，这一步具体的操作过程与步骤S2类似，在此不再赘述；
S4.3对获得的P波矢量波场和S波矢量波场进行各向异性介质弹性波逆时偏移，由震源端波场和检波器端波场，根据矢量成像条件把对应的波模互相关生成PP、PS、SP、SS偏移图像，并生成目标函数的梯度，梯度的计算公式为：
g＝LTr。
具体的，目标函数的梯度的计算公式为：
其中，g＝(gPP，gPS，gSP，gSS)T是目标函数的梯度，uP和φP分别是解耦后的震源端、检波器端P波矢量波场，uS和φS分别是解耦后的震源端、检波器端S波矢量波场；Tmax是地震数据最大接收时间。
S5计算步长和共轭梯度方向，结合梯度得到共轭梯度算法，采用共轭梯度算法对震源端波场进行更新迭代，重新合成多分量反射地震数据，计算目标函数，并判断目标函数是否满足收敛条件，直至目标函数满足收敛条件。即重复步骤S2-S5直至获得满足收敛条件。
步长α的计算公式为：
共轭梯度方向c(i+1)的计算公式为：
c(i+1)＝Pg(i+1)+βc(i)
其中，c(i)是第(i)次迭代的共轭梯度，L是各向异性弹性波逆时反偏移算子，P是预条件算子，β是系数。
共轭梯度算法为：
m(i+1)＝m(i)+αc(i+11，
其中，m(i)是第(i)次迭代的多分量弹性波图像模型，α是步长。
S6对各向异性介质弹性波高精度成像剖面进行绘图，并实现高分辨率图像输出，例如输出至电脑显示器，并在电脑磁盘中存储，但本实施例中进行图像显示的设备不止限于电脑，也可以是显示屏，手机，平板电脑等等。
实施例二
基于相同的发明构思，本实施例通过一个具体实例对本发明中各向异性介质弹性波高精度成像方法，进行说明，其具体步骤如下：
S1输入地层密度、垂向纵波横波速度、Thomsen各向异性参数的偏移背景模型，即各向异性参数模型和多分量弹性波图像(PP、PS、SP、SS)的初始值m(0)，如图4所示，模型离散为901×376个网格点，网格点间距为10米，该实验在地表均匀布设了181个炮和901个检波器，跑间距为50米，检波器间距为10米，实验的震源类型采用垂向力源，震源函数为主频15Hz的雷克子波，最大记录时间为6秒，时间采样间隔为1毫秒，观测数据的其中一个炮集记录如图5所示。通过交错网格有限差分，根据各向异性参数模型，对各向异性介质弹性波场正向传播数值进行模拟，求解各向异性弹性波动方程，得到震源端波场。该震源端波场如图6(a)和图6(b)所示。图6(a)是水平分量质点位移波场，图6(b)是垂直分量质点位移波场，震源端波场中纵波和横波耦合在一起，需要采用P/S波场矢量分解。
S2基于各向异性亥姆霍兹分解方法，对震源端波场进行矢量分解，获得的P波矢量波场和S波矢量波场，其结果如图6(c)、图6(d)、图6(e)和图6(f)所示，图6(c)和图6(d)分别是水平分量和垂直分量的纵波质点位移波场，得到图6(e)和图6(f)分别是水平分量和垂直分量的横波质点位移波场。获得的P波和S波矢量波场与对应的多分量弹性波图像相乘，得到波场-反射系数图像乘积，将波场-反射系数图像乘积分解成四种波模式，得到反射点处的PP波场，PS波场，SP波场和SS波场，将四种波模式结合作为二次震源进行正向传播数值模拟，在检波点处抽取由反偏移算子合成的反射波，合成多分量反射地震数据。
S3从电脑磁盘读入SEGY格式多分量实际地震数据炮集记录dobs，计算多分量反射地震数据与多分量实际地震数据dobs的目标函数，若目标函数满足收敛条件输出高分辨率图像，若目标函数不满足收敛条件，则进入下一步；
可选的，本实施例中的目标函数为：
其中，m(i)是当前第(i)次迭代的多分量弹性波图像；L是各向异性弹性波逆时反偏移算子，dobs是多炮多分量地震数据矢量；如果目标函数小于设定的阈值，则输出高分辨率图像，如果目标函数大于等于设定的阈值，则进入下一步。
S4计算多分量反射地震数据与多分量实际地震数据的残差r(i)＝dobs-Lm(i)，根据残差获得检波器端波场，根据震源端波场和检波器端波场获得目标函数的梯度；
弹性波逆时偏移的具体方法，如图3所示，为：
S4.1将残差作为伴随震源，通过波场反向传播模拟，即通过交错网格有限差分方法求解如下伴随波动方程得到检波器端波场φ，其结果如图7(a)和图7(b)所示，图7(a)是水平分量质点位移波场图，图7(b)是垂直分量质点位移波场图，震源端波场中纵波和横波耦合在一起，需要采用P/S波场矢量分解。
S4.2通过弹性波场矢量分解，将检波器端波场φ分解成P波矢量波场φP和S波矢量波场φS，其结果如图7(c)、图7(b)、图7(e)和图7(f)所示，图7(c)是解耦后的水平分量P波场图，图7(d)是解耦后的垂直分量P波场图，图7(e)是解耦后的水平分量S波场图，图7(f)为解耦后的垂直分量S波场图。
S4.3由震源端波场和检波器端波场，根据矢量成像条件把对应的波模互相关生成PP、PS、SP、SS偏移图像，并生成目标函数的梯度。
如图8所示，图8是现有技术中仅通过矢量成像条件把对应的波模互相关生成的VTI各向异性弹性波逆时偏移成像图；图8(a)是PP图，图8(b)是PS图，图8(c)是SP图，图8(d)是SS图；从图8中可知，通过现有技术中成像方法可以获得主要的反射界面，然而偏移图像的空间分辨率相对较低、存在低频噪声和波模式串扰，偏移图像振幅不均衡，不能准确地反映反射系数。
S5计算步长和共轭梯度方向，结合梯度得到共轭梯度算法，采用共轭梯度算法更新多分量弹性波图像模型，重新合成多分量反射地震数据，计算目标函数，并判断目标函数是否满足收敛条件，直至目标函数满足收敛条件。即重复步骤S2-S5直至获得满足收敛条件。
如图9所示，图9是基于本实施例中基于弹性波矢量分解的各向异性弹性波最小二乘逆时偏移生成的成像图；与图8对应，图9(a)是PP图，图9(b)是PS图，图9(c)是SP图，图9(d)是SS图。从图9中可知，通过本实施例中方法可以获得更加清晰的主要的反射界面，与图8中方案相比，本实施例中方案具有更高的空间分辨率，并通过过拟合数据压制了逆时偏移低频噪声。此外，通过最小二乘反演提高了图像振幅响应，深部图像的振幅得到了增强，高陡的盐丘侧面边界及盐下区域都得到了较好的成像效果。
如图10所示，图10展示了各向异性弹性波最小二乘逆时偏移目标函数收敛曲线，可以看出本实施例中方法基于弹性波场矢量分解的各向异性弹性波最小二乘逆时偏移具有良好的收敛性并能拟合数据。
S6对各向异性介质弹性波高精度成像剖面进行绘图，并实现高分辨率图像输出，例如输出至电脑显示器，并在电脑磁盘中存储，但本实施例中进行图像显示的设备不止限于电脑，也可以是显示屏，手机，平板电脑等等。
实施例三
基于相同的发明构思，本实施例公开了一种各向异性介质弹性波高精度成像系统，包括：
波场正向传播模拟模块，用于通过交错网格有限差分，根据各向异性参数模型，求解各向异性弹性波动方程，得到震源端波场；
弹性波场矢量分解模块，用于对所述震源端波场进行矢量分解，获得P波矢量波场和S波矢量波场；
各向异性弹性波逆时反偏移模块，用于通过震源端的所述P波矢量波场、S波矢量波场和输入的多分量弹性波图像初始模型，进行各向异性弹性波逆时反偏移，合成多分量反射地震数据；
目标函数判断模块，用于计算所述多分量反射地震数据与多分量实际地震数据的目标函数，若所述目标函数满足收敛条件输出高分辨率图像，若所述目标函数不满足收敛条件，则进入下一模块；
目标函数梯度计算模块，用于计算所述多分量反射地震数据与多分量实际地震数据的残差，将所述残差进行各向异性弹性波逆时偏移，获得目标函数的梯度；
共轭梯度法偏移图像更新模块，用于计算步长和共轭梯度方向，结合所述梯度得到共轭梯度算法，采用所述共轭梯度算法更新多分量弹性波图像模型，重新合成多分量反射地震数据，计算目标函数，并判断目标函数是否满足收敛条件，直至目标函数满足收敛条件。
基于上述技术方案，本发明具有以下技术效果：
1、本发明能够利用多分量地震勘探数据得到各向异性介质高精度多分量弹性波PP、PS、SP、SS图像剖面，为复杂油气藏提供更完备的成像信息，能更好地进行油气储层描述和地震解释；
2、本发明中方案把反演理论引入各向异性介质弹性波成像中，基于弹性波场矢量分解的各向异性弹性波最小二乘逆时偏移方法；相比传统的各向异性弹性波逆时偏移，本发明提高了图像的空间分辨率，减少了低频逆时偏移噪声，压制了波模式串扰噪声，显著提高了成像质量；
3、本发明的基于弹性波场矢量分解的各向异性弹性波最小二乘逆时偏移方法，基于最小二乘反演理论，在地震数据存在噪声或采样缺失时，能够嵌入数据权重函数，对质量差的数据赋予低的权重，减少对图像的影响。
4、本发明中方案为复杂油气藏精细勘探和开发提供支撑，助力对复杂油气藏的开发和利用，实现增储上产。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求保护范围之内。上述内容仅为本申请的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应以权利要求的保护范围为准。
Claims (10)
1.一种各向异性介质弹性波高精度成像方法，其特征在于，包括以下步骤：
通过交错网格有限差分，根据各向异性参数模型，求解各向异性弹性波动方程，得到震源端波场；
对所述震源端波场进行矢量分解，获得P波矢量波场和S波矢量波场；
利用震源端的所述P波矢量波场、S波矢量波场和输入的多分量弹性波图像初始模型，进行各向异性弹性波逆时反偏移，合成多分量反射地震数据；
计算所述多分量反射地震数据与多分量实际地震数据的目标函数，若所述目标函数满足收敛条件输出高分辨率图像，若所述目标函数不满足收敛条件，则进入下一步；
计算所述多分量反射地震数据与多分量实际地震数据的残差，将所述残差进行各向异性弹性波逆时偏移，获得目标函数的梯度；
计算步长和共轭梯度方向，结合所述梯度得到共轭梯度算法，采用所述共轭梯度算法更新多分量弹性波图像模型，重新合成多分量反射地震数据，计算目标函数，并判断目标函数是否满足收敛条件，直至目标函数满足收敛条件。
2.如权利要求1所述的各向异性介质弹性波高精度成像方法，其特征在于，所述目标函数为：
其中，m(i)是当前第(i)次迭代的多分量弹性波图像；L是各向异性弹性波逆时反偏移算子，dobs是多炮多分量地震数据矢量；如果目标函数小于设定的阈值，则输出高分辨率图像，如果目标函数大于等于设定的阈值，则进入下一步。
3.如权利要求1所述的各向异性介质弹性波高精度成像方法，其特征在于，所述震源端波场通过计算以下各向异性弹性波动方程获得：
其中，u＝(ux,uz)T是震源端波场，f＝(fx,fz)T是矢量震源，ρ是密度，Cij表示刚度矩阵的元素。
4.如权利要求3所述的各向异性介质弹性波高精度成像方法，其特征在于，对所述震源端波场进行矢量分解的方法为：通过快速傅里叶变换把弹性波场时间切片变换到波数-空间混合域，并求解各向异性矢量泊松方程，得到辅助矢量波场w，其计算公式为：
P波矢量波场uP(x,z)的计算公式为：
S波矢量波场uS(x,z)的计算公式为：
其中，是各向异性空间梯度算子。
5.如权利要求1所述的各向异性介质弹性波高精度成像方法，其特征在于，由各向异性弹性波逆时反偏移获得所述多分量反射地震数据的方法为：
将震源端波场进行矢量分解后获得的P波矢量波场和S波矢量波场与对应的波模式反射系数图像相乘，得到波场-反射系数图像乘积；
将所述波场-反射系数图像乘积分解成四种波模式，得到反射点处的PP波场，PS波场，SP波场和SS波场；
将所述四种波模式结合作为二次震源，对所述二次震源进行基于交错网格有限差分的正向传播数值模拟，在检波点处抽取由反偏移算子合成的反射波，合成多分量反射地震数据。
6.如权利要求1所述的各向异性介质弹性波高精度成像方法，其特征在于，所述各向异性弹性波逆时偏移的具体方法为：
将所述残差作为伴随震源，通过波场反向传播模拟，得到检波器端波场；
通过弹性波场矢量分解，将所述检波器端波场分解成P波矢量波场和S波矢量波场；
由所述震源端波场和检波器端波场，根据矢量成像条件把对应的波模互相关生成PP、PS、SP、SS偏移图像，并生成目标函数的梯度。
7.如权利要求6所述的各向异性介质弹性波高精度成像方法，其特征在于，所述检波器端波通过计算以下异性弹性波动方程获得：
其中，φ＝(φx,φz)T是检波器端波场，r＝(rx,rz)T是伴随震源，ρ是密度，Cij表示刚度矩阵的元素。
8.如权利要求6所述的各向异性介质弹性波高精度成像方法，其特征在于，所述目标函数的梯度的公式为：
其中，g＝(gPP,gPS,gSP,gSS)T是目标函数的梯度，uP和φP分别是解耦后的震源端、检波器端P波矢量波场，uS和φS分别是解耦后的震源端、检波器端S波矢量波场；Tmax是地震数据的最大接收时间。
9.如权利要求8所述的各向异性介质弹性波高精度成像方法，其特征在于，所述步长α的计算公式为：
共轭梯度方向c(i+1)的计算公式为：
c(i+1)＝Pg(i+1)+βc(i)
其中，L是各向异性弹性波逆时反偏移算子，P是预条件算子，β是系数。
10.一种各向异性介质弹性波高精度成像系统，其特征在于，包括：
波场正向传播模拟模块，用于通过交错网格有限差分，根据各向异性参数模型，求解各向异性弹性波动方程，得到震源端波场；
弹性波场矢量分解模块，用于对所述震源端波场进行矢量分解，获得P波矢量波场和S波矢量波场；
各向异性弹性波逆时反偏移模块，用于通过震源端的所述P波矢量波场、S波矢量波场和输入的多分量弹性波图像初始模型，进行各向异性弹性波逆时反偏移，合成多分量反射地震数据；
目标函数判断模块，用于计算所述多分量反射地震数据与多分量实际地震数据的目标函数，若所述目标函数满足收敛条件输出高分辨率图像，若所述目标函数不满足收敛条件，则进入下一模块；
目标函数梯度计算模块，用于计算所述多分量反射地震数据与多分量实际地震数据的残差，将所述残差进行各向异性弹性波逆时偏移，获得目标函数的梯度；
共轭梯度法偏移图像更新模块，用于计算步长和共轭梯度方向，结合所述梯度得到共轭梯度算法，采用所述共轭梯度算法更新多分量弹性波图像模型，重新合成多分量反射地震数据，计算目标函数，并判断目标函数是否满足收敛条件，直至目标函数满足收敛条件。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310105843.7A CN116413801B (zh) | 2023-02-13 | 2023-02-13 | 一种各向异性介质弹性波高精度成像方法和系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310105843.7A CN116413801B (zh) | 2023-02-13 | 2023-02-13 | 一种各向异性介质弹性波高精度成像方法和系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116413801A CN116413801A (zh) | 2023-07-11 |
| CN116413801B true CN116413801B (zh) | 2023-09-29 |
Family
ID=87052236
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310105843.7A Active CN116413801B (zh) | 2023-02-13 | 2023-02-13 | 一种各向异性介质弹性波高精度成像方法和系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116413801B (zh) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105974470A (zh) * | 2016-07-04 | 2016-09-28 | 中国石油大学(华东) | 一种多分量地震资料最小二乘逆时偏移成像方法及系统 |
| CN107589443A (zh) * | 2017-08-16 | 2018-01-16 | 东北石油大学 | 基于弹性波最小二乘逆时偏移成像的方法及系统 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101549388B1 (ko) * | 2014-10-17 | 2015-09-02 | 한국지질자원연구원 | 탄성파 다성분 자료에 대한 중합전 ｅｇｓ 구조보정 방법 |
| US11733413B2 (en) * | 2021-04-30 | 2023-08-22 | Saudi Arabian Oil Company | Method and system for super resolution least-squares reverse time migration |
-
2023
- 2023-02-13 CN CN202310105843.7A patent/CN116413801B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105974470A (zh) * | 2016-07-04 | 2016-09-28 | 中国石油大学(华东) | 一种多分量地震资料最小二乘逆时偏移成像方法及系统 |
| CN107589443A (zh) * | 2017-08-16 | 2018-01-16 | 东北石油大学 | 基于弹性波最小二乘逆时偏移成像的方法及系统 |
Non-Patent Citations (1)
| Title |
|---|
| 何兵寿 ; 张会星 ; .多分量波场的矢量法叠前深度偏移技术.石油地球物理勘探.2006,(第04期),全文. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116413801A (zh) | 2023-07-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11092707B2 (en) | Determining a component of a wave field | |
| Ren et al. | A hierarchical elastic full-waveform inversion scheme based on wavefield separation and the multistep-length approach | |
| US8830788B2 (en) | Sensitivity kernal-based migration velocity analysis in 3D anisotropic media | |
| EP2497043B1 (en) | Seismic imaging systems and methods employing a 3d reverse time migration with tilted transverse isotropy | |
| Yu et al. | Acoustic-elastic coupled equation for ocean bottom seismic data elastic reverse time migration | |
| US20120016592A1 (en) | Image domain signal to noise estimate with borehole data | |
| Wu et al. | Microseismic source locations with deconvolution migration | |
| CN113740901B (zh) | 基于复杂起伏地表的陆上地震数据全波形反演方法及装置 | |
| Sun et al. | Elastic full-waveform inversion with extrapolated low-frequency data | |
| Gao et al. | An efficient vector elastic reverse time migration method in the hybrid time and frequency domain for anisotropic media | |
| Li et al. | Waveform inversion of seismic first arrivals acquired on irregular surface | |
| Ma et al. | Seismic full‐waveform inversion of the crust‐mantle structure beneath China and adjacent regions | |
| Li et al. | Continent‐continent collision between the south and north China plates revealed by seismic refraction and reflection at the southern segment of the Tanlu fault zone | |
| CN116413801B (zh) | 一种各向异性介质弹性波高精度成像方法和系统 | |
| CN114814944B (zh) | 基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法 | |
| Ren et al. | Pre-stack elastic reverse time migration in tunnels based on cylindrical coordinates | |
| CN115373022A (zh) | 一种基于振幅相位校正的弹性波场Helmholtz分解方法 | |
| MX2011003852A (es) | Atributos de procesamiento de imagen con inversion de tiempo. | |
| CN111208568A (zh) | 一种时间域多尺度全波形反演方法及系统 | |
| Yu et al. | Analytic acoustic–elastic coupled equations and their application in vector-wave-based imaging of ocean bottom 4C data | |
| Roberts et al. | Investigation into the use of 2D elastic waveform inversion from look‐ahead walk‐away VSP surveys | |
| Chi et al. | Source-independent amplitude-semblance full-waveform inversion using a hybrid time-and frequency-domain approach | |
| Zheng et al. | Improved numerical solution of anisotropic poroelastic wave equation in microseismicity: Graphic process unit acceleration and moment tensor implementation | |
| Plessix et al. | Frequency-domain finite-difference migration with only few frequencies? | |
| Poletto et al. | Seismic acquisition and processing of onshore dual fields by a reciprocal experiment |
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 |