CN110858005B - 一种基于基追踪横向多道约束的各向异性参数反演方法 - Google Patents

一种基于基追踪横向多道约束的各向异性参数反演方法 Download PDF

Info

Publication number
CN110858005B
CN110858005B CN201810975237.XA CN201810975237A CN110858005B CN 110858005 B CN110858005 B CN 110858005B CN 201810975237 A CN201810975237 A CN 201810975237A CN 110858005 B CN110858005 B CN 110858005B
Authority
CN
China
Prior art keywords
parameter
inversion
inverted
matrix
transverse
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
Application number
CN201810975237.XA
Other languages
English (en)
Other versions
CN110858005A (zh
Inventor
刘婵娟
胡华锋
许凯
梁志强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China Petroleum and Chemical Corp, Sinopec Geophysical Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201810975237.XA priority Critical patent/CN110858005B/zh
Publication of CN110858005A publication Critical patent/CN110858005A/zh
Application granted granted Critical
Publication of CN110858005B publication Critical patent/CN110858005B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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/40—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging
    • G01V1/44—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging using generators and receivers in the same well
    • G01V1/48—Processing data
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/40—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging
    • G01V1/44—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging using generators and receivers in the same well
    • G01V1/48—Processing data
    • G01V1/50—Analysing data
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/61—Analysis by combining or comparing a seismic data set with other data
    • G01V2210/616—Data from specific type of measurement
    • G01V2210/6169—Data from specific type of measurement using well-logging
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/622—Velocity, density or impedance
    • G01V2210/6226—Impedance
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/626—Physical property of subsurface with anisotropy

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Remote Sensing (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种基于基追踪横向多道约束的各向异性参数反演方法，包括：处理地震数据，获取方位角度叠加道集；构建岩石物理模型，获取裂缝参数测井数据；对所述测井数据进行去相关处理；构建反演目标函数；基于所述反演目标函数，获得待反演参数反射系数序列。本发明的基于基追踪横向多道约束的各向异性参数反演方法，通过纳入方位角获得方位角度叠加道集，方位角度叠加道集的信道比高，且包括多重地震信息，进而较为准确地提供储层中的裂缝信息，通过对测井数据进行去相关处理，提高了反演结果的准确度，通过在目标函数中加入相关约束项，构建反演目标函数，提高了反演结果稳定性和横向连续性，以此降低未知参数估测的不确定性。

Description

一种基于基追踪横向多道约束的各向异性参数反演方法
技术领域
本发明属于地球物理反演领域，具体涉及一种基于基追踪横向多道约束的各向异性参数反演方法。
背景技术
随着勘探目标转向复杂型油气藏，非常规油气勘探开发已成为研究重点。在页岩气的勘探中，有效识别裂缝发育程度是至关重要的，而裂缝型储层具有各向异性的特点，因此，准确获得能够反映裂缝发育程度的裂缝各向异性参数对地下裂缝预测具有重要意义。
随着宽方位地震数据采集和处理技术的快速发展，利用宽方位叠前地震数据预测地下裂缝已成为主要研究方法。目前基于方位叠前地震数据的各向异性参数反演方法主要利用地震波在裂缝型储层传播过程中反射系数随入射角和方位角的变化情况来进行研究，其中理论推导主要来源于Ruger在1998年提出的基于等效横向各向同性介质的反射系数近似公式，随后有相关学者对该方程进行了一系列简化，但在各向异性参数反演中仍存在可靠性和稳定性的问题。另外裂缝型储层的各向异性特征对地震波振幅产生的影响，使得在反演过程中出现横向不连续的情况。对此特别需要一种合理的各向异性参数反演方法，能提高反演的可靠性、稳定性和横向连续性，从而有助于开展利用地震信息测量地下各向异性强弱的方法研究，进而较为准确地提供储层中的裂缝信息，指导裂缝型储层的勘探和开发。
发明内容
本发明的目的是提出一种提高反演准确性、可靠性、稳定性和横向连续性的基于基追踪横向多道约束的各向异性参数反演方法。
为了实现上述目的，本发明提供一种基于基追踪横向多道约束的各向异性参数反演方法，包括：处理地震数据，获取方位角度叠加道集；构建裂缝岩石物理模型，获取裂缝参数测井数据；对所述测井数据进行去相关处理；构建反演目标函数；基于所述反演目标函数，获得待反演参数反射系数序列。
优选的，所述对所述测井数据进行去相关处理包括：生成所述测井数据的协方差矩阵；对所述协方差矩阵进行奇异值分解，获得去相关矩阵；基于所述去相关矩阵，对公式d＝Gr进行变换，获得变换后的去相关公式，其中，d为所述方位角度叠加道集的地震数据，G为子波矩阵，r为待反演参数。
优选的，所述变换后的去相关公式为：d＝G×r，其中，G＝G×V，r＝V-1×r，d为所述方位角度叠加道集的地震数据，G为去相关的子波矩阵，G为子波矩阵，V为去相关矩阵，r为去相关的待反演参数，r为待反演参数。
优选的，所述构建反演目标函数包括：获得所述方位角度叠加道集的似然函数；基于所述似然函数，获得目标函数；在所述目标函数中加入相关约束项，构建所述反演目标函数。
优选的，所述基于所述似然函数，获得目标函数包括：结合所述似然函数与先验分布，获得后验概率分布；对所述后验概率分布进行最大化，获得所述目标函数。
优选的，所述相关约束项包括：基于井数据平滑阻抗模型约束、L1范数约束和横向多道约束。
优选的，所述反演目标函数为：
Figure BDA0001777245200000031
其中，d为所述方位角度叠加道集的地震数据，G为子波矩阵，r为待反演参数，Jsmooth(r)为基于井数据平滑阻抗模型约束，r'为前一道待反演参数反演结果，λ、β为权系数。
优选的，通过基追踪算法获得所述待反演参数反射系数序列。
优选的，通过对所述反演目标函数求取最小值获得所述待反演参数反射系数序列。
优选的，根据所述待反演参数反射系数序列，获得最终待反演参数，其中，所述最终待反演参数包括纵波阻抗、横波阻抗、裂缝岩石物理参数。
本发明的有益效果在于：本发明的基于基追踪横向多道约束的各向异性参数反演方法，通过对地震数据进行处理获得方位角度叠加道集，通过构建岩石物理模型，获取裂缝参数测井数据，对测井数据之间的协方差矩阵进行去相关处理，构建反演目标函数，从而获得待反演参数反射系数序列，通过纳入方位角获得方位角度叠加道集，方位角度叠加道集的信道比高，且包括多重地震信息，进而较为准确地提供储层中的裂缝信息，指导裂缝型储层的勘探和开发。
本发明的基于基追踪横向多道约束的各向异性参数反演方法，通过对测井数据进行去相关处理，提高了反演结果的准确度，为地下裂缝预测提供较为准确的信息。
本发明的基于基追踪横向多道约束的各向异性参数反演方法，通过在目标函数中加入相关约束项，构建反演目标函数，再对反演目标函数求取获得待反演参数反射系数序列，提高了反演结果稳定性和横向连续性，以此降低未知参数估测的不确定性，为地下裂缝预测提供合理可靠且有力的证据，解决了裂缝储层各向异性参数反演中存在的不稳定、横向不连续等问题，对于指导裂缝有效识别具有显著意义。
本发明具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施例中将是显而易见的，或者将在并入本文中的附图和随后的具体实施例中进行详细陈述，这些附图和具体实施例共同用于解释本发明的特定原理。
附图说明
通过结合附图对本发明示例性实施方式进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施方式中，相同的参考标号通常代表相同部件。
图1示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的流程图。
图2示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的方位角小角度叠加数据。
图3示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的方位角中角度叠加数据。
图4示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的方位角大角度叠加数据。
图5示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的纵波阻抗反演结果。
图6示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的横波阻抗反演结果。
图7示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的裂缝岩石物理参数Γx反演结果。
图8示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的裂缝岩石物理参数Γy反演结果。
图9a示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的纵波阻抗IP反演结果对比。
图9b示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的横波阻抗IS反演结果对比。
图9c示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的裂缝岩石物理参数Γx反演结果对比。
图9d示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的裂缝岩石物理参数Γy反演结果对比。
图10a示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的纵波阻抗Ip横向多道对比。
图10b示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的横波阻抗IS横向多道对比。
图10c示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的裂缝岩石物理参数Γx横向多道对比。
图10d示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的裂缝岩石物理参数Γy横向多道对比。
具体实施方式
下面将更详细地描述本发明的优选实施方式。虽然以下描述了本发明的优选实施方式，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。
根据本发明的基于基追踪横向多道约束的各向异性参数反演方法，包括：处理地震数据，获取方位角度叠加道集；构建裂缝岩石物理模型，获取裂缝参数测井数据；对测井数据进行去相关处理；构建反演目标函数；基于反演目标函数，获得待反演参数反射系数序列。
具体地，对地震数据进行处理获得方位角度叠加道集，通过构建岩石物理模型，获取裂缝参数测井数据，测井数据包括：纵波阻抗、横波阻抗、裂缝岩石物理参数，再对测井数据之间进行去相关处理，构建反演目标函数，从而获得待反演参数反射系数序列。
根据示例性的实施方式基于基追踪横向多道约束的各向异性参数反演方法，通过纳入方位角获得方位角度叠加道集，方位角度叠加道集的信道比高，且包括多重地震信息，进而较为准确地提供储层中的裂缝信息，指导裂缝型储层的勘探和开发。
作为优选方案，对测井数据进行去相关处理包括：生成测井数据的协方差矩阵；对协方差矩阵进行奇异值分解，获得去相关矩阵；基于去相关矩阵，对公式d＝Gr进行变换，获得变换后的去相关公式，其中，d为方位角度叠加道集的地震数据，G为子波矩阵，r为待反演参数。
作为优选方案，变换后的去相关公式为：d＝G×r，其中，G＝G×V，r＝V-1×r，d为方位角度叠加道集的地震数据，G为去相关的子波矩阵，G为子波矩阵，V为去相关矩阵，r为去相关的待反演参数，r为待反演参数。
具体的，由于地层的纵波阻抗、横波阻抗和裂缝岩石物理参数之间是统计相关的，需要应用测井数据参数之间的协方差矩阵对测井数据进行去相关处理。首先将测井数据的协方差矩阵表示成：
Figure BDA0001777245200000061
其中：
Figure BDA0001777245200000062
为纵波阻抗反射系数的方差，
Figure BDA0001777245200000063
为纵波阻抗反射系数与横波阻抗反射系数之间的协方差，
Figure BDA0001777245200000064
为纵波阻抗反射系数与岩石物理参数之间的协方差，σIs 2为横波阻抗反射系数的方差，
Figure BDA0001777245200000065
为横波阻抗反射系数与岩石物理参数之间的协方差，
Figure BDA0001777245200000071
为岩石物理参数之间的协方差，
Figure BDA0001777245200000072
为岩石物理参数的方差。
一般来说，由于这四个变量间相互不独立，协方差矩阵Cm的非对角线上的元素不为零，对该协方差矩阵进行奇异值分解：
Figure BDA0001777245200000073
其中，
Figure BDA0001777245200000074
为纵波阻抗的方差，
Figure BDA0001777245200000075
为横波阻抗的方差，
Figure BDA0001777245200000076
为岩石物理参数的方差。
假设
Figure BDA0001777245200000077
在假定反射系数不变且独立的情况下将其扩展为N个时间采样，则去相关矩阵V：
Figure BDA0001777245200000078
对公式d＝Gr进行变换：
Figure BDA0001777245200000079
得：
d＝G×r (5)
其中，d为方位角度叠加道集的地震数据，G为子波矩阵，G为去相关的子波矩阵，V为去相关矩阵，r为待反演参数，r为去相关的待反演参数。
由公式(5)可以得出，去相关处理后的d方位角度叠加道集的地震数据与原来的d方位角度叠加道集的地震数据相同，因此，去相关处理对方位角度叠加道集的地震数据没有影响。
所以，经过变换后的四参数协方差矩阵Cm'等于：
Figure BDA0001777245200000081
可以看出，Cm'非对角线上的元素都为零，表明经过变换后的四个参数之间是相互独立的。
根据示例性的实施方式的基于基追踪横向多道约束的各向异性参数反演方法，通过对测井数据进行去相关处理，提高了反演结果的准确度，为地下裂缝预测提供较为准确的信息。
作为优选方案，构建反演目标函数包括：获得所述方位角度叠加道集的似然函数；基于似然函数，获得目标函数；在目标函数中加入相关约束项，构建反演目标函数。
具体的，假设地震资料背景噪声服从高斯分布，则可得到方位角度叠加道集的似然函数：
Figure BDA0001777245200000082
其中，d为所述方位角度叠加道集的地震数据，G为子波矩阵，r为待反演参数，σn 2是噪音的方差。
基于似然函数获得目标函数，同时将去相关处理后获得的参数代入目标函数中。
通常通过对目标函数进行反演从而获得可靠的待反演参数，但是仅靠地震数据与反演结果间的误差最小二乘拟合不能得到精确反演解，在此基础上在目标函数中加入了相关约束项来增加反演稳定性以及横向连续性，从而减少反演结果的多解性，进而得到有意义的地层反射系数序列，最终获得可靠的待反演参数。
根据示例性的实施方式的基于基追踪横向多道约束的各向异性参数反演方法，通过在目标函数中加入相关约束项，构建反演目标函数，再对反演目标函数求取获得待反演参数反射系数序列，提高了反演结果稳定性和横向连续性，以此降低未知参数估测的不确定性，为地下裂缝预测提供合理可靠且有力的证据，解决了裂缝储层各向异性参数反演中存在的不稳定、横向不连续等问题，对于指导裂缝有效识别具有显著意义。
作为优选方案，基于似然函数，获得目标函数包括：结合似然函数与先验分布，获得后验概率分布；对后验概率分布进行最大化，获得目标函数。
具体的，将似然函数与以前的数据经验结合起来得到后验概率分布，将后验概率分布最大化，得到目标函数：
J(r)＝(d-Gr)T(d-Gr) (8)
其中，d为所述方位角度叠加道集的地震数据，G为子波矩阵，r为待反演参数。
作为优选方案，相关约束项包括：基于井数据平滑阻抗模型约束、L1范数约束和横向多道约束。
具体的，基于井数据平滑阻抗模型约束是在有井资料时，通过数学变换把测井信息加入到反演方程中进一步约束反射系数，使反演问题变得更加适定，更有意义，同时，也可以地提高了反演的精度和分辨率。
地下地层反射系数为稀疏序列，通过加入L1范数约束使反射系数稀疏化，反演结果更稳定。
在地层平稳沉积中，地层横向基本连续，表现在地震反演中即道与道之间具有相关性，通过加入道与道之间的约束，即横向多道约束，来提高横向连续性。
作为优选方案，反演目标函数为：
Figure BDA0001777245200000101
其中，
Figure BDA0001777245200000102
d为方位角度叠加道集的地震数据，G为子波矩阵，r为待反演参数，Jsmooth(r)为基于井数据平滑阻抗模型约束，r'为前一道待反演参数反演结果，λ、β为权系数，αp为纵波阻抗反射系数约束对应的系数项，αs为横波阻抗反射系数约束对应的系数项，
Figure BDA0001777245200000103
为裂缝岩石物理参数反射系数约束对应的系数项，Cp为纵波阻抗对应的积分矩阵，Cs为横波阻抗对应的积分矩阵，
Figure BDA0001777245200000104
为裂缝岩石物理参数对应的积分矩阵，ξP为纵波阻抗初始模型，ξs为横波阻抗初始模型，
Figure BDA0001777245200000105
Figure BDA0001777245200000106
为裂缝岩石物理参数初始模型，E为单位矩阵，ξ'为前一道待反演参数反演结果。
具体的，Jsmooth(r)为基于井数据平滑阻抗模型约束，λ||r||1为L1范数约束，
Figure BDA0001777245200000111
为横向多道约束。
在地震数据与反演结果间的误差最小二乘拟合下，使用公式r＝V-1×r，其中，V为去相关矩阵，r为去相关的待反演参数，r为待反演参数，把测井数据的纵波阻抗、横波阻抗、裂缝岩石物理参数这四个待反演参数乘以V-1，获得去相关处理后的待反演参数加入到反演目标函数中进一步约束反射系数。
以纵波阻抗为例，一般有：
Figure BDA0001777245200000112
其中，Ip(t)为纵波阻抗，ΔIp(t)为上下地层纵波阻抗差，t为某个采样点。
对反射系数积分有：
Figure BDA0001777245200000113
其中，t＝[t0,...,tM-1]T，Ip(t0)为初始纵波波阻抗值，rp(η)为反射系数。
写成矩阵形式，即：
ξp＝Cprp (12)
其中，rp为反射系数，
Figure BDA0001777245200000114
且
Figure BDA0001777245200000115
为积分算子矩阵，其离散形式可表示为
Figure BDA0001777245200000116
所以基于井数据平滑阻抗模型约束项可表示为：
Figure BDA0001777245200000121
其中，αp为纵波阻抗反射系数约束对应的系数项，αs为横波阻抗反射系数约束对应的系数项，
Figure BDA0001777245200000122
为裂缝岩石物理参数反射系数约束对应的系数项，Cp为纵波阻抗对应的积分矩阵，Cs为横波阻抗对应的积分矩阵，
Figure BDA0001777245200000123
为裂缝岩石物理参数对应的积分矩阵，ξP为纵波阻抗初始模型，ξs为横波阻抗初始模型，
Figure BDA0001777245200000124
为裂缝岩石物理参数初始模型。
反演目标函数可表示为：
Figure BDA0001777245200000125
其中，
Figure BDA0001777245200000126
d为方位角度叠加道集的地震数据，G为子波矩阵，r为待反演参数，Jsmooth(r)为基于井数据平滑阻抗模型约束，r'为前一道待反演参数反演结果，λ、β为权系数，αp为纵波阻抗反射系数约束对应的系数项，αs为横波阻抗反射系数约束对应的系数项，
Figure BDA0001777245200000127
为裂缝岩石物理参数反射系数约束对应的系数项，Cp为纵波阻抗对应的积分矩阵，Cs为横波阻抗对应的积分矩阵，
Figure BDA0001777245200000128
为裂缝岩石物理参数对应的积分矩阵，ξP为纵波阻抗初始模型，ξs为横波阻抗初始模型，
Figure BDA0001777245200000129
Figure BDA0001777245200000131
为裂缝岩石物理参数初始模型，E为单位矩阵，ξ'为前一道待反演参数反演结果。。
作为优选方案，通过基追踪算法获得待反演参数反射系数序列。
作为优选方案，通过对反演目标函数求取最小值获得待反演参数反射系数序列。
作为优选方案，根据待反演参数反射系数序列，获得最终待反演参数，其中，最终待反演参数包括纵波阻抗、横波阻抗、裂缝岩石物理参数。
具体的，对得到的目标函数求取最小值得到待反演参数反射系数序列，进而推导得到最终待反演参数，其中包括纵波阻抗、横波阻抗、裂缝岩石物理参数。
实施例
图1示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的流程图。
如图1所示，基于基追踪横向多道约束的各向异性参数反演方法，包括：
S102：处理地震数据，获取方位角度叠加道集；
S104：构建裂缝岩石物理模型，获取裂缝参数测井数据；
根据实际工期及经验公式构建裂缝岩石物理模型，基于裂缝岩石物理模型获取裂缝参数与弹性参数之间的关系，得到裂缝参数的响应，即获得了裂缝参数测井数据，其中，测井数据包括：纵波阻抗、横波阻抗、裂缝岩石物理参数。
S106：对测井数据进行去相关处理；
其中，对测井数据进行去相关处理包括：生成测井数据的协方差矩阵；对协方差矩阵进行奇异值分解，获得去相关矩阵；基于去相关矩阵，对公式d＝Gr进行变换，获得变换后的去相关公式，其中，d为方位角度叠加道集的地震数据，G为子波矩阵，r为待反演参数。
其中，变换后的去相关公式为：d＝G×r，其中，G＝G×V，r＝V-1×r，d为方位角度叠加道集的地震数据，G为去相关的子波矩阵，G为子波矩阵，V为去相关矩阵，r为去相关的待反演参数，r为待反演参数。
具体的，由于地层的纵波阻抗、横波阻抗和裂缝岩石物理参数之间是统计相关的，需要应用测井数据参数之间的协方差矩阵对测井数据进行去相关处理。首先将测井数据的协方差矩阵表示成：
Figure BDA0001777245200000141
其中：
Figure BDA0001777245200000142
为纵波阻抗反射系数的方差，
Figure BDA0001777245200000143
为纵波阻抗反射系数与横波阻抗反射系数之间的协方差，
Figure BDA0001777245200000144
为纵波阻抗反射系数与岩石物理参数之间的协方差，σIs 2为横波阻抗反射系数的方差，
Figure BDA0001777245200000145
为横波阻抗反射系数与岩石物理参数之间的协方差，
Figure BDA0001777245200000146
为岩石物理参数之间的协方差，
Figure BDA0001777245200000147
为岩石物理参数的方差。
一般来说，由于这四个变量间相互不独立，协方差矩阵Cm的非对角线上的元素不为零，对该协方差矩阵进行奇异值分解：
Figure BDA0001777245200000148
其中，
Figure BDA0001777245200000149
为纵波阻抗的方差，
Figure BDA00017772452000001410
为横波阻抗的方差，
Figure BDA00017772452000001411
为岩石物理参数的方差。
假设
Figure BDA0001777245200000151
在假定反射系数不变且独立的情况下将其扩展为N个时间采样，则去相关矩阵V：
Figure BDA0001777245200000152
对公式d＝Gr进行变换：
Figure BDA0001777245200000153
得：
d＝G×r (5)
其中，d为方位角度叠加道集的地震数据，G为子波矩阵，G为去相关的子波矩阵，V为去相关矩阵，r为待反演参数，r为去相关的待反演参数。
由公式(5)可以得出，去相关处理后的d方位角度叠加道集的地震数据与原来的d方位角度叠加道集的地震数据相同，因此，去相关处理对方位角度叠加道集的地震数据没有影响。
所以，经过变换后的四参数协方差矩阵Cm'等于：
Figure BDA0001777245200000154
可以看出，Cm'非对角线上的元素都为零，表明经过变换后的四个参数之间是相互独立的。
S108：构建反演目标函数；
S108构建反演目标函数包括子步骤S1081-S1084：
S1081：假设地震资料背景噪声服从高斯分布，则可得到方位角度叠加道集的似然函数：
Figure BDA0001777245200000161
其中，d为方位角度叠加道集的地震数据，G为子波矩阵，r为待反演参数，σn 2是噪音的方差。
S1082：将似然函数与先验分布结合起来得到后验概率分布，将后验概率分布最大化，得到目标函数：
J(r)＝(d-Gr)T(d-Gr) (8)
S1083：将去相关处理后获得的参数代入目标函数。
S1084：在目标函数中加入相关约束项，构建反演目标函数。
一般来说，仅靠地震数据与反演结果间的误差最小二乘拟合不能得到精确反演解，在此基础上本发明在目标函数中加入了相关约束项来增加反演稳定性以及横向连续性，从而减少反演结果的多解性，进而得到有意义的地层反射系数序列，最终获得可靠的待反演参数。这里的相关约束项包括基于井数据平滑阻抗模型约束、L1范数约束和横向多道约束，最终获得的反演目标函数为：
Figure BDA0001777245200000162
其中Jsmooth(r)为基于井数据平滑阻抗模型约束项，λ||r||1为L1范数约束，
Figure BDA0001777245200000163
为横向多道约束，r'为前一道待反演参数反演结果，λ、β为权系数。
具体的，使用公式r＝V-1×r，其中，V为去相关矩阵，r为去相关的待反演参数，r为待反演参数，把测井数据的纵波阻抗、横波阻抗、裂缝岩石物理参数这四个待反演参数乘以V-1，获得去相关处理后的待反演参数加入到反演目标函数中进一步约束反射系数。
以纵波阻抗为例，一般有：
Figure BDA0001777245200000171
对反射系数积分有：
Figure BDA0001777245200000172
其中，t＝[t0,...,tM-1]T，Ip(t0)为初始纵波波阻抗值。写成矩阵形式，即：
ξp＝Cprp (12)
其中，
Figure BDA0001777245200000173
且
Figure BDA0001777245200000174
为积分算子矩阵，其离散形式可表示为
Figure BDA0001777245200000175
所以基于井数据平滑阻抗模型约束项可表示为：
Figure BDA0001777245200000176
反演目标函数可表示为：
Figure BDA0001777245200000177
其中，
Figure BDA0001777245200000181
d为方位角度叠加道集的地震数据，G为子波矩阵，r为待反演参数，Jsmooth(r)为基于井数据平滑阻抗模型约束，r'为前一道待反演参数反演结果，λ、β为权系数，αp为纵波阻抗反射系数约束对应的系数项，αs为横波阻抗反射系数约束对应的系数项，
Figure BDA0001777245200000183
Figure BDA0001777245200000184
为裂缝岩石物理参数反射系数约束对应的系数项，Cp为纵波阻抗对应的积分矩阵，Cs为横波阻抗对应的积分矩阵，
Figure BDA0001777245200000182
为裂缝岩石物理参数对应的积分矩阵，ξP为纵波阻抗初始模型，ξs为横波阻抗初始模型，
Figure BDA0001777245200000185
为裂缝岩石物理参数初始模型，E为单位矩阵，ξ'为前一道待反演参数反演结果。
S110：基于反演目标函数，获得待反演参数反射系数序列。
通过基追踪算法对反演目标函数求取最小值获得待反演参数反射系数序列。
最后的步骤，根据待反演参数反射系数序列，获得最终待反演参数，其中，最终待反演参数包括纵波阻抗、横波阻抗、裂缝岩石物理参数。
该基于基追踪横向多道约束的各向异性参数反演方法的工作过程如下：对地震数据进行处理得到方位角度叠加道集，根据工区构建岩石物理模型进而计算获取裂缝参数测井数据，由于地层的纵波阻抗、横波阻抗和裂缝岩石物理参数之间是统计相关的，应用它们之间的协方差矩阵对测井数据进行去相关处理，去相关处理时首先生成待反演的弹性参数及裂缝岩石物理参数之间的协方差矩阵，对该协方差矩阵进行奇异值分解，获得去相关矩阵V，再对公式d＝Gr进行变换，得到变换后的去相关公式d＝G×r，可知去相关处理对方位角度叠加道集的地震数据没有影响；基于方位角度叠加道集生成似然函数，将似然函数结合以前的数据经验，获得后验概率分布，对后验概率分布进行最大化，从而获得目标函数，然后将测井数据参数乘以V-1，获得去相关处理后的参数，将去相关后的参数代入目标函数，再将基于井数据平滑阻抗模型约束、L1范数约束和横向多道约束的参数也加入目标函数，最后获得反演目标函数，通过基追踪算法对反演目标函数求取最小值获得待反演参数反射系数序列，获得最终待反演参数，其中，最终待反演参数包括纵波阻抗、横波阻抗、裂缝岩石物理参数。
图2、图3和图4分别示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的方位角小角度叠加数据、方位角中角度叠加数据、方位角大角度叠加数据。在图2-4中，横坐标表示地震道数，纵坐标表示采样点数。
如图2、图3、图4所示，对地震数据信息通过纳入方位角获得的方位角度叠加道集，该方位角度叠加道集反应的是真实的地震数据。
图5、图6分别示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的纵波阻抗IP反演结果、横波阻抗IS反演结果；图7、图8分别示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的裂缝岩石物理参数Γx反演结果、裂缝岩石物理参数Γy反演结果。在图5-8中，横坐标表示地震道数，纵坐标表示采样点数。
如图5、图6、图7和图8所示，通过基于基追踪横向多道约束的各向异性参数反演方法得到的纵波阻抗IP反演结果、横波阻抗IS反演结果、裂缝岩石物理参数Γx反演结果，以及裂缝岩石物理参数Γy反演结果，由图可以看出纵波阻抗IP反演结果、横波阻抗IS反演结果、裂缝岩石物理参数Γx反演结果，以及裂缝岩石物理参数Γy反演结果都与地震数据吻合，基本能反映裂缝发育程度。
图9a、图9b、图9c和图9d分别示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的纵波阻抗IP、横波阻抗IS、裂缝岩石物理参数Γx、裂缝岩石物理参数Γy的反演结果对比。
图10a、图10b、图10c和图10d分别示出了根据本发明的一个实施例的基于基追踪横向多道约束的各向异性参数反演方法的纵波阻抗Ip、横波阻抗IS、裂缝岩石物理参数Γx、裂缝岩石物理参数Γy的横向多道对比。
如9a所示，通过基于基追踪横向多道约束的各向异性参数反演方法抽取的纵波阻抗IP反演结果对比，黑色点线为反演结果，灰色直线为模型数据，横坐标的数值是纵波阻抗(kg/m3)·(m/s)，由图可以看出纵波阻抗Ip反演结果与模型数据比较吻合。如9b所示，通过基于基追踪横向多道约束的各向异性参数反演方法抽取的横波阻抗IS反演结果对比，黑色点线为反演结果，灰色直线为模型数据，横坐标的数值是横波阻抗(kg/m3)·(m/s)，由图可以看出横波阻抗IS反演结果稍不理想，这也是下一步需要改进的地方。如9c所示，通过基于基追踪横向多道约束的各向异性参数反演方法抽取的裂缝岩石物理参数Γx反演结果对比，黑色点线为反演结果，灰色直线为模型数据，由图可以看出裂缝岩石物理参数Γx反演结果与模型数据比较吻合。如9d所示，通过基于基追踪横向多道约束的各向异性参数反演方法抽取的裂缝岩石物理参数Γy反演结果对比，黑色点线为反演结果，灰色直线为模型数据，由图可以看出裂缝岩石物理参数Γy反演结果与模型数据比较吻合。在图9a，纵坐标表示采样点数，横坐标表示纵波阻抗(kg/m2*s)，在图9b，纵坐标表示采样点数，横坐标表示横波阻抗(kg/m2*s)，在图9c，纵坐标表示采样点数，横坐标表示岩石物理参数τx，在图9d，纵坐标表示采样点数，横坐标表示岩石物理参数τy。
图10a所示，通过基于基追踪横向多道约束的各向异性参数反演方法抽取的纵波阻抗Ip横向多道对比，黑色点线为反演结果，灰色直线为模型数据，横坐标的数值是纵波阻抗(kg/m3)·(m/s)，由图可以看出纵波阻抗IP横向上反演结果也吻合模型数据，横向连续性也得到了提高。图10b所示，通过基于基追踪横向多道约束的各向异性参数反演方法抽取的横波阻抗IS横向多道对比，黑色点线为反演结果，灰色直线为模型数据，横坐标的数值是横波阻抗(kg/m3)·(m/s)，由图可以看出横波阻抗IS横向上反演结果也吻合模型数据，横向连续性也得到了提高。如10c所示通过基于基追踪横向多道约束的各向异性参数反演方法抽取的裂缝岩石物理参数Γx横向多道对比，黑色点线为反演结果，灰色直线为模型数据，由图可以看出裂缝岩石物理参数Γx横向上反演结果也吻合模型数据，横向连续性也得到了提高。如图10d所示，通过基于基追踪横向多道约束的各向异性参数反演方法抽取的裂缝岩石物理参数Γy横向多道对比，黑色点线为反演结果，灰色直线为模型数据，由图可以看出裂缝岩石物理参数Γy横向上反演结果也吻合模型数据，横向连续性也得到了提高。在图10a，纵坐标表示纵波阻抗(kg/m2*s)，横坐标表示采样点数，在图10b，纵坐标表示横波阻抗(kg/m2*s)，横坐标表示采样点数，在图10c，纵坐标表示岩石物理参数τx，横坐标表示采样点数，在图10d，纵坐标表示岩石物理参数τy，横坐标表示采样点数。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。

Claims (7)

1.一种基于基追踪横向多道约束的各向异性参数反演方法，其特征在于，包括：
处理地震数据，获取方位角度叠加道集；
构建裂缝岩石物理模型，获取裂缝参数测井数据；
对所述测井数据进行去相关处理；
构建反演目标函数；
基于所述反演目标函数，获得待反演参数反射系数序列；
所述构建反演目标函数包括：
获得所述方位角度叠加道集的似然函数；
基于所述似然函数，获得目标函数；
在所述目标函数中加入相关约束项，构建所述反演目标函数；
所述反演目标函数为：
Figure FDA0003227771580000011
其中，d为方位角度叠加道集的地震数据，G为子波矩阵，r为待反演参数，Jsmooth(r)为基于井数据平滑阻抗模型约束，r'为前一道待反演参数反演结果，λ、β为权系数，λ||r||1为L1范数约束，
Figure FDA0003227771580000012
为横向多道约束；
其中，基于井数据平滑阻抗模型约束项可表示为：
Figure FDA0003227771580000013
αp为纵波阻抗反射系数约束对应的系数项，αs为横波阻抗反射系数约束对应的系数项，
Figure FDA0003227771580000014
为裂缝岩石物理参数反射系数约束对应的系数项，Cp为纵波阻抗对应的积分矩阵，Cs为横波阻抗对应的积分矩阵，
Figure FDA0003227771580000015
Figure FDA0003227771580000016
为裂缝岩石物理参数对应的积分矩阵，ξP为纵波阻抗初始模型，ξs为横波阻抗初始模型，
Figure FDA0003227771580000021
为裂缝岩石物理参数初始模型。
2.根据权利要求1所述的各向异性参数反演方法，其特征在于，所述对所述测井数据进行去相关处理包括：
生成所述测井数据的协方差矩阵；
对所述协方差矩阵进行奇异值分解，获取去相关矩阵；
基于所述去相关矩阵，对公式d＝Gr进行变换，获得变换后的去相关公式，
其中，d为所述方位角度叠加道集的地震数据，G为子波矩阵，r为待反演参数。
3.根据权利要求2所述的各向异性参数反演方法，其特征在于，所述变换后的去相关公式为：
d＝G×r
其中，G＝G×V，r＝V-1×r，d为所述方位角度叠加道集的地震数据，G为去相关的子波矩阵，G为子波矩阵，V为去相关矩阵，r为去相关的待反演参数，r为待反演参数。
4.根据权利要求1所述的各向异性参数反演方法，其特征在于，所述基于所述似然函数，获得目标函数包括：
结合所述似然函数与先验分布，获得后验概率分布；
对所述后验概率分布进行最大化，获得所述目标函数。
5.根据权利要求1所述的各向异性参数反演方法，其特征在于，通过基追踪算法获得所述待反演参数反射系数序列。
6.根据权利要求5所述的各向异性参数反演方法，其特征在于，通过对所述反演目标函数求取最小值获得所述待反演参数反射系数序列。
7.根据权利要求1所述的各向异性参数反演方法，其特征在于，还包括：
根据所述待反演参数反射系数序列，获得最终待反演参数，其中，所述最终待反演参数包括纵波阻抗、横波阻抗、裂缝岩石物理参数。
CN201810975237.XA 2018-08-24 2018-08-24 一种基于基追踪横向多道约束的各向异性参数反演方法 Active CN110858005B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810975237.XA CN110858005B (zh) 2018-08-24 2018-08-24 一种基于基追踪横向多道约束的各向异性参数反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810975237.XA CN110858005B (zh) 2018-08-24 2018-08-24 一种基于基追踪横向多道约束的各向异性参数反演方法

Publications (2)

Publication Number Publication Date
CN110858005A CN110858005A (zh) 2020-03-03
CN110858005B true CN110858005B (zh) 2021-11-05

Family

ID=69636380

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810975237.XA Active CN110858005B (zh) 2018-08-24 2018-08-24 一种基于基追踪横向多道约束的各向异性参数反演方法

Country Status (1)

Country Link
CN (1) CN110858005B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114076980B (zh) * 2020-08-17 2024-04-02 中国石油化工股份有限公司 一种针对薄层刻画的方法及系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102788994A (zh) * 2012-07-12 2012-11-21 恒泰艾普石油天然气技术服务股份有限公司 一种储层裂缝的确定方法
CN104005760A (zh) * 2014-04-16 2014-08-27 孙赞东 基于方位各向异性弹性阻抗的裂缝检测方法
CN104007467A (zh) * 2014-04-16 2014-08-27 孙赞东 一种基于混合范数正则化的叠前三参数反演实现的储层与流体预测方法
CN106199695A (zh) * 2016-06-29 2016-12-07 中国石油化工股份有限公司 基于空变目标函数的叠前三参数反演方法
CN107783188A (zh) * 2016-08-31 2018-03-09 中国石油化工股份有限公司 一种针对岩石力学参数的叠前反演方法及装置

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7519477B2 (en) * 2006-10-20 2009-04-14 Bhp Billiton Innovation Pty Ltd. Method for determining impedance coefficients of a seismic trace
GB2555365B (en) * 2015-07-28 2021-07-28 Geoquest Systems Bv Seismic constrained discrete fracture network

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102788994A (zh) * 2012-07-12 2012-11-21 恒泰艾普石油天然气技术服务股份有限公司 一种储层裂缝的确定方法
CN104005760A (zh) * 2014-04-16 2014-08-27 孙赞东 基于方位各向异性弹性阻抗的裂缝检测方法
CN104007467A (zh) * 2014-04-16 2014-08-27 孙赞东 一种基于混合范数正则化的叠前三参数反演实现的储层与流体预测方法
CN106199695A (zh) * 2016-06-29 2016-12-07 中国石油化工股份有限公司 基于空变目标函数的叠前三参数反演方法
CN107783188A (zh) * 2016-08-31 2018-03-09 中国石油化工股份有限公司 一种针对岩石力学参数的叠前反演方法及装置

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
Multitrace impedance inversion with lateral constraints;Haitham Hamid et al;《GEOPHYSICS》;20151231;第80卷(第6期);第M102-M103页 *
Stochastic inversion of elastic impedance based on simulated annealing genetic algorithm;Chanjuan Liu et al;《2016 SEG International Exposition and 86th Annual Meeting》;20161231;第2896-2900页 *
叠前三参数同步反演方法及其应用;杨培杰等;《石油学报》;20090331;第30卷(第2期);第233-234页 *
基于贝叶斯理论的孔隙流体模量叠前AVA反演;李超等;《石油物探》;20150731;第54卷(第4期);第467-476页 *
基于贝叶斯线性AVAZ的TTI介质裂缝参数反演;葛子建等;《地球物理学报》;20180731;第61卷(第7期);第3010-3013页 *
非均匀正交各向异性特征参数地震散射反演方法研究;潘新朋等;《地球物理学报》;20180131;第61卷(第1期);第267-281页 *

Also Published As

Publication number Publication date
CN110858005A (zh) 2020-03-03

Similar Documents

Publication Publication Date Title
US8103453B2 (en) Method of seismic data interpolation by projection on convex sets
Křížová et al. Resolvability of isotropic component in regional seismic moment tensor inversion
CA2683618C (en) Inverse-vector method for smoothing dips and azimuths
CN108459350A (zh) 一种深度域地震子波提取与地震记录合成的一体化方法
CN110456417B (zh) 一种地震数据多次波压制方法
CN103954992B (zh) 一种反褶积方法及装置
Wang et al. Seismic velocity inversion transformer
CN110895348B (zh) 一种地震弹性阻抗低频信息提取方法、系统及存储介质
WO2014084751A1 (en) A method for processing acoustic waveforms
CN103869362B (zh) 体曲率获取方法和设备
Aleardi et al. Characterisation of shallow marine sediments using high‐resolution velocity analysis and genetic‐algorithm‐driven 1D elastic full‐waveform inversion
MacBeth et al. Methods of measurement for 4D seismic post‐stack time shifts
Dettmer et al. Joint time/frequency-domain inversion of reflection data for seabed geoacoustic profiles and uncertainties
CN110858005B (zh) 一种基于基追踪横向多道约束的各向异性参数反演方法
Aeron et al. Broadband dispersion extraction using simultaneous sparse penalization
CN112946752B (zh) 基于裂缝密度模型预测裂缝概率体的方法
Wang et al. Deep learning‐based H‐κ method (HkNet) for estimating crustal thickness and Vp/Vs ratio from receiver functions
Jiang et al. Full waveform inversion based on inversion network reparameterized velocity
Guigné et al. Acoustic zoom high-resolution seismic beamforming for imaging specular and non-specular energy of deep oil and gas bearing geological formations
CN111273346B (zh) 去除沉积背景的方法、装置、计算机设备及可读存储介质
Wu et al. An unsupervised inversion method for seismic brittleness parameters driven by the physical equation
CN111983668B (zh) 一种获取地震参数估计的方法及系统
CN116068644A (zh) 一种利用生成对抗网络提升地震数据分辨率和降噪的方法
CN114325832A (zh) 一种裂缝参数和弹性参数同步反演方法及系统
CN111077573A (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