CN114063156A - 基于稀疏表征的地震各向异性梯度反演方法及系统 - Google Patents

基于稀疏表征的地震各向异性梯度反演方法及系统 Download PDF

Info

Publication number
CN114063156A
CN114063156A CN202010748575.7A CN202010748575A CN114063156A CN 114063156 A CN114063156 A CN 114063156A CN 202010748575 A CN202010748575 A CN 202010748575A CN 114063156 A CN114063156 A CN 114063156A
Authority
CN
China
Prior art keywords
inversion
sparse
intermediate variable
gradient
anisotropic
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.)
Granted
Application number
CN202010748575.7A
Other languages
English (en)
Other versions
CN114063156B (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.)
Petrochina Co Ltd
Original Assignee
Petrochina Co Ltd
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 Petrochina Co Ltd filed Critical Petrochina Co Ltd
Priority to CN202010748575.7A priority Critical patent/CN114063156B/zh
Priority claimed from CN202010748575.7A external-priority patent/CN114063156B/zh
Publication of CN114063156A publication Critical patent/CN114063156A/zh
Application granted granted Critical
Publication of CN114063156B publication Critical patent/CN114063156B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/306Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/40Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/15Correlation function computation including computation of convolution operations
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/61Analysis by combining or comparing a seismic data set with other data
    • G01V2210/616Data from specific type of measurement
    • G01V2210/6161Seismic or acoustic, e.g. land or sea measurements
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/61Analysis by combining or comparing a seismic data set with other data
    • G01V2210/616Data from specific type of measurement
    • G01V2210/6169Data from specific type of measurement using well-logging
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/62Physical property of subsurface
    • G01V2210/626Physical property of subsurface with anisotropy

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Remote Sensing (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geophysics (AREA)
  • Acoustics & Sound (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Geology (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种基于稀疏表征的地震各向异性梯度反演方法及系统,方法包括:根据测井数据得到中间变量及各向异性梯度;根据中间变量建立训练集,通过K‑SVD方法从训练集中学习中间变量对应的稀疏字典;对中间变量进行插值处理,建立各向异性介质中各向异性梯度反演的初始模型;根据各向异性梯度理论和稀疏表征理论建立目标函数,其中,假设目标函数在同一工区中且地层在横向上具有连续性,将中间变量对应的稀疏字典上的稀疏表达作为约束条件;将初始模型作为初始值,利用目标函数进行反演计算,得到中间变量的反演结果,并根据各向异性梯度与中间变量的关系得到各向异性梯度的反演结果,该反演结果趋向于分辨率更高的测井数据,准确性更高。

Description

基于稀疏表征的地震各向异性梯度反演方法及系统
技术领域
本发明属于地震数据解释领域,尤指一种基于稀疏表征的地震各向异性梯度反演方法及系统。
背景技术
地震各向异性梯度反演是利用地表观测的叠前方位角道集地震数据,以已知地质规律和钻井、测井资料为约束,对地下岩层裂缝及各向异性特征的空间结构和进行成像的过程。各向异性梯度反演为油气勘探中的储层预测提供了重要的依据,利用地表观测的宽方位地震数据可以更好的反演各向异性梯度。
各向异性梯度反演框架同叠前各向同性反演一致,是一种基于模型的叠前反演方法,是预测地层裂缝储层的重要手段。该反演方法从人为给定的初始地质模型出发,利用合成地震记录与实际地震数据之间的残差计算出模型的修正量,以此来不断迭代更新初始模型,使模型正演合成记录与实际地震记录达到最佳吻合,最终得到的模型数据便是反演结果。
各向异性梯度反演时有四个反演参数,参数多且参数对地震数据影响微弱等原因导致反演结果存在极大地不稳定性。通常解决反演稳定性最直接方法就是加入先验信息约束,但是传统的先验约束信息都需要假设模型具有特定的结构(稀疏、光滑等)或具有一个唯一确定的数学表达形式(高斯或者柯西先验分布)。比如,增强反演结果光滑性最流行的正则化方法是Tikhonov型正则化,该类方法以最小化模型参数的低阶空间导数的l2范数作为先验约束;以及,2012年,Anagaw等基于全变差(TV)方法设计出能生成块状结构的弹性参数的反演方法,而这种方法是以最小化模型参数的梯度的l1范数为目标。这些先验信息通常只是一种单纯的数学抽象,只能适用于特定的工区,对地下信息的刻画并不完备,如果将这种先验信息用于各向异性反演梯度必然导致反演结果存在误差。
由于传统的正则化方法没有考虑到测井数据中内含的地层特征,针对这一缺陷,基于稀疏表征的反演方法因此产生。稀疏表征技术是将信号在预先获得的一组基向量中进行分解,使得信号可以尽可能地由少量几个基向量线性表示。如今,稀疏表征理论在图像超分辨率重建,信号去噪,等领域中都取得成功的应用。
稀疏表征技术包括字典学习和稀疏重构。1996年,Olshausen等在《Nature》上发表的著名的Sparsenet字典学习方法奠定了字典学习的理论基础。其主要思想是从自然图像库中抽取大量的小图像块作为训练集,采用交替优化方法不断的学习训练更新字典,以获得信号的稀疏表征。2006年,Aharon提出了K-SVD方法,可很好地用于字典学习中;该方法使用正交匹配追踪技术,逐个考虑原子信号对数据拟合的贡献并逐个更新原子信号,目前,该方法已经成为字典学习应用最广泛的方法。稀疏重构就是将信号在预先学好的字典中进行稀疏分解,使得信号可以尽可能地由少量原子线性表示。稀疏分解实质是求一个l0范数的问题,求解方法主要有匹配追踪,正交匹配追踪等。
字典学习最早用于主成分分析方法中,2006年,Elad成功的将字典学习应用于图像去噪。2012年,Tang等引入了学习型超完备字典,将字典学习与总变差(TV)最小化结合起来用于地震数据去噪,并证明非均匀原子的字典更适合于地震数据去噪。2015年,Zhu改进了以往的字典方法,引入基于参数的字典学习方法用于地震数据去噪,在合成地震记录中有很好的效果。2016年,Chen将字典模型的分析方法和学习方法结合起来提出双稀疏字典学习方法用于地震数据去噪。在2017年,Nazari等通过多任务策略将二维的字典学习模型扩展到三维并用于三维地震数据去噪。同样地,2018年,Liu等将字典学习用于地震数据的稀疏表征,进而作为一个正则项用于地震数据去噪。前述在地震勘探中的字典学习方法均是用于数据去噪。2018年,Li首次将字典学习方法用于地震的全波形反演,提出基于多类字典学习的正则化框架,并在数值模拟中有较好的成果。同样地,在2018年,She等首次将字典学习方法运用于叠前AVO弹性参数反演,并提出将稀疏表征作为正则化项这一新方法,不管是在数值模拟还是在实际数据上均有很好的效果,这是反演历史的一个突破;但是,该方法只能解决各向同性的弹性参数反演,不能对各向异性介质中的地震数据进行反演。
综上来看,亟需一种可以解决上述问题的能够对地震各向异性梯度进行反演的技术方案。
发明内容
为解决上述问题,本发明提出了一种基于稀疏表征的地震各向异性梯度反演方法及系统。该方法及系统可以保证各向异性梯度反演中的稳定性,能够更准确的发现地下的各向异性情况,使得反演结果更趋向于分辨率更高的测井数据,进而保证纵向分辨率的提高,使得反演结果更准确。
在本发明实施例的第一方面,提出了一种基于稀疏表征的地震各向异性梯度反演方法,该方法包括:
采集测井数据;
根据所述测井数据得到中间变量及各向异性梯度;
根据所述中间变量建立训练集,通过K-SVD方法从所述训练集中学习中间变量对应的稀疏字典;
对所述中间变量进行插值处理,建立各向异性介质中各向异性梯度反演的初始模型;
根据各向异性梯度理论和稀疏表征理论建立目标函数,其中,假设所述目标函数在同一工区中且地层在横向上具有连续性,将所述中间变量对应的稀疏字典上的稀疏表达作为约束条件;
将所述初始模型作为初始值,利用所述目标函数进行反演计算,得到中间变量的反演结果,并根据所述各向异性梯度与中间变量的关系得到各向异性梯度的反演结果。
在本发明实施例的第二方面,提出了一种基于稀疏表征的地震各向异性梯度反演系统,该系统包括:
数据采集模块,用于采集测井数据;
参数计算模块,用于根据所述测井数据得到中间变量及各向异性梯度;
字典训练模块,用于根据所述中间变量建立训练集,通过K-SVD方法从所述训练集中学习中间变量对应的稀疏字典;
插值处理模块,用于对所述中间变量进行插值处理,建立各向异性介质中各向异性梯度反演的初始模型;
目标函数建立模块,用于根据各向异性梯度理论和稀疏表征理论建立目标函数,其中,假设所述目标函数在同一工区中且地层在横向上具有连续性,将所述中间变量对应的稀疏字典上的稀疏表达作为约束条件;
反演计算模块,用于将所述初始模型作为初始值,利用所述目标函数进行反演计算,得到中间变量的反演结果,并根据所述各向异性梯度与中间变量的关系得到各向异性梯度的反演结果。
在本发明实施例的第三方面,提出了一种计算机设备,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行所述计算机程序时实现基于稀疏表征的地震各向异性梯度反演方法。
在本发明实施例的第四方面,提出了一种计算机可读存储介质,所述计算机可读存储介质存储有计算机程序,所述计算机程序被处理器执行时实现基于稀疏表征的地震各向异性梯度反演方法。
本发明提出的基于稀疏表征的地震各向异性梯度反演方法及系统,能够更准确的发现地下的各向异性情况,根据实际测井数据获得中间变量,基于K-SVD方法从实际数据中学习得中间变量参数各自的稀疏字典,并将反演的各参数在对应的稀疏字典上的稀疏表达作为约束条件,从而进行各参数的反演,所得到的反演结果更趋向于分辨率更高的测井数据,实现纵向分辨率的提高且反演结果更准确,为后续的储层预测提供有力的数据支持。
附图说明
为了更清楚地说明本申请实施例技术方案,下面将对实施例描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图是本申请的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其它的附图。
图1是本发明一实施例的基于稀疏表征的地震各向异性梯度反演方法流程示意图。
图2是本发明一实施例的反演计算的详细流程示意图。
图3A及图3B分别是本发明一具体实施例的利用测井数据训练得到的中间变量A和B的前40个原子的示意图。
图4A及图4B分别是本发明一具体实施例的中间变量A经过反演的结果和经过Tikh正则化的反演结果示意图。
图5A及图5B分别是本发明一具体实施例的中间变量B经过反演的结果和经过Tikh正则化的反演结果示意图。
图6A及图6B分别是本发明一具体实施例的各向异性梯度Bani经过反演的结果和经过Tikh正则化的反演结果示意图。
图7是本发明一实施例的基于稀疏表征的地震各向异性梯度反演系统架构示意图。
图8是本发明一实施例的计算机设备结构示意图。
具体实施方式
下面将参考若干示例性实施方式来描述本发明的原理和精神。应当理解,给出这些实施方式仅仅是为了使本领域技术人员能够更好地理解进而实现本发明,而并非以任何方式限制本发明的范围。相反,提供这些实施方式是为了使本公开更加透彻和完整,并且能够将本公开的范围完整地传达给本领域的技术人员。
本领域技术人员知道,本发明的实施方式可以实现为一种系统、装置、设备、方法或计算机程序产品。因此,本公开可以具体实现为以下形式,即:完全的硬件、完全的软件(包括固件、驻留软件、微代码等),或者硬件和软件结合的形式。
根据本发明的实施方式,提出了一种基于稀疏表征的地震各向异性梯度反演方法及系统。
下面参考本发明的若干代表性实施方式,详细阐释本发明的原理和精神。
图1是本发明一实施例的基于稀疏表征的地震各向异性梯度反演方法流程示意图。如图1所示,该方法包括:
步骤S101,采集测井数据;
步骤S102,根据所述测井数据得到中间变量及各向异性梯度;
步骤S103,根据所述中间变量建立训练集,通过K-SVD方法从所述训练集中学习中间变量对应的稀疏字典;
步骤S104,对所述中间变量进行插值处理,建立各向异性介质中各向异性梯度反演的初始模型;
步骤S105,根据各向异性梯度理论和稀疏表征理论建立目标函数,其中,假设所述目标函数在同一工区中且地层在横向上具有连续性,将所述中间变量对应的稀疏字典上的稀疏表达作为约束条件;
步骤S106,将所述初始模型作为初始值,利用所述目标函数进行反演计算,得到中间变量的反演结果,并根据所述各向异性梯度与中间变量的关系得到各向异性梯度的反演结果。
在一实施例中,步骤S102,根据所述测井数据得到中间变量A、B、C、D及各向异性梯度Bani,其中,
Figure BDA0002609256360000061
Figure BDA0002609256360000062
C=Bani cos 2φiso
D=Bani sin 2φiso
Figure BDA0002609256360000063
Figure BDA0002609256360000064
其中,A、B、C、D为中间变量;Bani为各向异性梯度;Biso为各向同性梯度;φiso为入射角;ΔZ为上下界面的纵波阻抗之差;
Figure BDA0002609256360000065
为上下界面的纵波阻抗的均值;Δδ、Δγ分别为上下界面的各项异性参数之差;
Figure BDA0002609256360000066
为上下界面的纵波速度平均值;
Figure BDA0002609256360000067
为上下界面的横波速度平均值;Δα为上下界面的纵波速度之差;ΔG为上下界面的横波剪切模量之差;
Figure BDA0002609256360000068
为上下界面的横波剪切模量平均值。
在一实施例中,步骤S103包括:
根据所述中间变量A、B、C、D建立训练集,其中训练集中包含多个训练样本;
通过K-SVD方法从所述训练集中学习中间变量对应的稀疏字典,供反演计算使用。
在一实施例中,参考图2,为本发明一实施例的反演计算的详细流程示意图。如图2所示,步骤S106包括:
步骤S1061,将初始模型中的中间变量进行分块处理,得到各参数的N个长度为M的训练样本集;
步骤S1062,设置目标函数中的正则化参数以及最大迭代次数;
步骤S1063,以所述中间变量对应的稀疏字典作为中间变量对应的初始字典,根据所述初始字典确定对应的初始化的稀疏系数矩阵;
步骤S1064,根据正则化参数、最大迭代次数及目标函数,基于正则化方法反演得到中间变量的结果;
步骤S1065,将稀疏字典中的原子固定,采用OMP方法对训练样本集进行稀疏编码,更新初始化的稀疏系数矩阵;
步骤S1066,根据更新后的稀疏系数矩阵及中间变量的结果,重构出对应反演层段的变量值,得到中间变量的反演结果,并根据所述各向异性梯度与中间变量的关系得到各向异性梯度的反演结果。
为了对上述基于稀疏表征的地震各向异性梯度反演方法进行更为清楚的解释,下面结合一个具体的实施例来进行说明。
本发明提出的基于稀疏表征的地震各向异性梯度反演方法主要包含两部分:字典学习和稀疏反演。
在字典学习阶段,先通过测井数据计算得到中间变量A、B、C、D,并且学习对应变量的字典,即由测井数据得到的中间变量分块得到训练数据,再运用K-SVD等方法学得稀疏字典,该字典可表述测井数据中蕴含的地层特征;其中,在She等2018年公开的文献中,对于利用测井数据进行字典学习的过程进行了详细的描述。
字典学习和稀疏表征在学术界的称谓是稀疏字典学习;稀疏表征理论最早是在研究信号处理应用中发展起来的。其基础是多尺度分析理论,在此基础上拓展,形成了相应的理论框架。字典学习的最终任务是从许多信号中推知原子信号以及各个信号的原子信号权重,它的实质是对大量数据信号的一种降维表示,它尝试找出隐藏在信号数据中每个信号所共有的特征。稀疏表征则是可以从原子信号中重构出原始信号,即用尽可能少的原子信号来表示原始信号,这种表示还能带来一个附加的好处,即计算速度快。近年来稀疏字典学习在图像处理领域中得到广泛应用。以下为本发明采稀疏字典学习的具体实现过程:
在根据测井数据训练字典时,先通过测井数据计算得到中间变量A、B、C、D。本发明需要将中间变量分成小块(小块是一个一维的数据),以作为原始信号,每个信号的长度为M,M的取值范围一般在40到120之间,可根据实际情况确定。
利用这些小块合并,作为训练样本集;设训练样本集为:
Figure BDA0002609256360000081
其中,R为Y中训练样本的个数。
字典学习过程的目标是找到一个字典
Figure BDA0002609256360000082
能以最稀疏的方式表达训练集Y,若记此时的稀疏系数集为:
Figure BDA0002609256360000083
则该过程可以被表示为如下能量最小化问题:
Figure BDA0002609256360000084
其中,每一个长度为L的列向量αi为样本yi在字典D下的稀疏表征;||·||0是l0-范数,用于计算某一个向量的含零元素的个数;K是表示稀疏程度的常量。在式(1)中,代价函数
Figure BDA0002609256360000085
确立了对训练集的分解的一致性,而约束项则限制了表达结果的稀疏性。
基于K-SVD的测井数据的初始字典学习方法的流程如下:
S1,对已知的每条测井曲线的数据分别进行分块处理,得到不同参数的R个长度为M的训练样本Y=[y1 y2…yM];
S2,设置正则化参数λ,最大迭代次数T;
S3,初始化字典D和稀疏系数矩阵α,设置当前迭代次数t=1;
S4,稀疏编码:假设D中的原子固定,用OMP方法对训练样本集Y进行稀疏编码,得到稀疏系数矩阵α。
S5,字典更新:对字典中的原子进行逐个更新。对于第k(k=1,2,…,K)个原子,按照下述步骤更新:
S51,dk为当前要进行更新的原子,记Ik={i|αi(k)≠0,1≤i≤N},αi(k)为αi中的第k个元素,则Ik表示所有样本中使用到第k列原子dk的索引集。
S52,用
Figure BDA0002609256360000086
表示去除字典中第k个原子后的字典,
Figure BDA0002609256360000087
表示去除第k行系数以后的系数矩阵,然后计算新的字典与系数矩阵与原信号之间的残差矩阵
Figure BDA0002609256360000088
Ek表示除去第k个原子后对所有样本造成的误差;
S53,根据Ik中的索引值选取Ek中相应的列向量,构成新的误差矩阵
Figure BDA0002609256360000089
并对
Figure BDA00026092563600000810
进行奇异值分解
Figure BDA00026092563600000811
S54,取矩阵U中的第一列,即最大特征值对应的特征向量,作为更新后的原子,并更新相应的系数。
S6,设置t'=t+1。检查t'是否大于T。如果t≤T,则返回步骤S4;如果t>T,则停止迭代,输出字典D。
构建AVAZ反演框架:
各向异性假设下双层HTI介质中纵波反射系数近似表达式;
Figure BDA0002609256360000091
其中,i、φ分别是入射角和方位角;α、β分别为纵波和横波速度;Z为纵波阻抗,Z=ρa;G为横波剪切模量,G=ρβ2;运算符Δ[·]为上下界面的参数之差,
Figure BDA0002609256360000092
为上下界面参数之均值,ρ为密度,ε、δ、γ为各向异性参数。
当入射角较小时,式(2)关于入射角的高阶项(第三项)可以被忽略,因此可以将式(2)简化为:
Figure BDA0002609256360000093
其中,
Figure BDA0002609256360000094
Figure BDA0002609256360000095
Figure BDA0002609256360000096
A为中间变量;Bani为各向异性梯度;Biso为各向同性梯度;ΔZ为上下界面的纵波阻抗之差;
Figure BDA0002609256360000097
为上下界面的纵波阻抗的均值;Δδ、Δγ分别为上下界面的各项异性参数之差;
Figure BDA0002609256360000098
为上下界面的纵波速度平均值;
Figure BDA0002609256360000099
为上下界面的横波速度平均值;Δα为上下界面的纵波速度之差;ΔG为上下界面的横波剪切模量之差;
Figure BDA00026092563600000910
为上下界面的横波剪切模量平均值。
基于式(3),对非线性方位反射系数公式线性化,可以得到:
Figure BDA0002609256360000101
其中,
Figure BDA0002609256360000102
C=Bani cos 2φiso (6)
D=Bani sin 2φiso
B、C、D为中间变量;Bani为各向异性梯度;Biso为各向同性梯度;φiso为入射角。
因此,可将各向异性梯度用下式计算:
Figure BDA0002609256360000103
所以,通过上述推导可以得到中间变量A、B、C、D,然后再由式(7)可以计算各向异性梯度Bani,用以表征裂缝发育情况。
单道地震记录d(t)可以通过子波w(t)与反射系数卷积而得:
Figure BDA0002609256360000104
建立数学模型如下:
d=Em (9)
其中,E表示由子波和权系数组成的核矩阵;m表示模型参数。
Figure BDA0002609256360000105
其中,N为地质模型采样点数,A1,…,AN表示中间变量A。
根据式(8),已知叠前方位纵波地震数据d,和矩阵E,估算地下参数m,此过程称为反演,通常,该过程目标函数如下:
Figure BDA0002609256360000106
其中,λ是正则化参数;
利用现有的正则化方法来约束反演;比如,Tikhonov正则化和全变差正则化。Tikhonov正则化通过最小化模型m的梯度的二范使得解具有光滑性:
Figure BDA0002609256360000107
全变差正则化通过假设模型m具有稀疏的梯度,使得解呈块状特征:
Figure BDA0002609256360000111
式(12)、(13)中T为一阶差分算子。
由于上述正则化约束存在各向异性参数反演中的稳定性差的问题,对此,本发明提出了一种能够根据数据本身构建自适应的先验模型的方法,以实现更各向异性参数更精确的反演,降低反演不稳定性。
由于地壳运动的尺度很大,故同一工区中,地层通常在横向上具有连续性。上述通过测井学得的字典使得每条测井数据都能在该字典上稀疏表达。由于同一工区中,地层通常在横向上具有连续性,故对于每个地震道的反演结果通常也能在该字典上稀疏表达。根据这条思路,将模型在字典上的稀疏表达作为约束条件来约束反演。
将式(10)重写如下:
Figure BDA0002609256360000112
其中,n是经过块处理后模型的块的数量,以中间向量A为例,
Figure BDA0002609256360000113
表示从上述稀疏字典学习过程学得的具有L个原子的字典,
Figure BDA0002609256360000114
表示训练样本的第i个小块,它是长度为M的一个列向量,相应的,每一个列向量
Figure BDA0002609256360000115
是通过字典DA
Figure BDA0002609256360000116
的稀疏表征结果。
同理,上述描述对中间变量B、C、D也是一样的。
此外,μ1、μ2、μ3、μ4表示四个合适的权重参数。Ri是一个4M×4N的矩阵,用于提取模型参数m的第i个小块。
分析目标函数,第一项表示反演数据和真实数据的误差,后面四项则为反演结果在字典上稀疏表达的整体误差。约束项是使得反演参数能在各自对应字典上进行稀疏表征。求解(11)式由于参数太多过于困难,四个参数必须得同时考虑。因此,本发明引入中间变量x,将目标函数写成式(15),约束项不变,式中β是一个新的正则化参数。
Figure BDA0002609256360000121
采用坐标下降方法,将求解式(15),求解过程分为以下两步:
S01,求解目标函数(16):
Figure BDA0002609256360000122
当k=0时,设置x0=m0(初始模型),该步骤是一个凸优化问题,求解方法众多,如梯度下降方法、拟牛顿方法等。
S02,求解目标函数(17):
Figure BDA0002609256360000123
式(17)等价于对每个参数求形如式(16)的目标函数:
Figure BDA0002609256360000131
其中,j=1,2,3,4,依次代表需要反演的四个参数。
经过上述过程,利用拟牛顿方法和OMP方法最小化目标函数,可以得到地层各向异性参数的反演结果。
在上述过程中,有多个参数需要设定,比如原子长度和正则化参数,这些参数的选取需要做质量控制,用不同的参数去反演,根据反演效果评估确定最优参数。
需要说明的是,尽管在上述实施例及附图中以特定顺序描述了本发明方法的操作,但是,这并非要求或者暗示必须按照该特定顺序来执行这些操作,或是必须执行全部所示的操作才能实现期望的结果。附加地或备选地,可以省略某些步骤,将多个步骤合并为一个步骤执行,和/或将一个步骤分解为多个步骤执行。
进一步的,为了对上述基于稀疏表征的地震各向异性梯度反演方法进行更为清楚的解释,下面结合一个具体的实施例来进行说明,然而值得注意的是该实施例仅是为了更好地说明本发明,并不构成对本发明不当的限定。
以中国西部某工区为例,采集该工区的实际数据,利用本发明的方法进行反演测试。
根据测井数据计算中间变量A、B、C、D,并分别训练中间变量共四个参数的稀疏字典,如图3A和图3B所示,分别为利用测井数据训练得到的中间变量A和B的前40个原子,这些原子代表了测井数据代表的不同地质特征。
利用本发明的基于稀疏表征的地震各向异性梯度反演方法,应用于该工区的一条剖面,反演结果如图4A至6B所示;图4A、图4B分别显示了中间变量A经过反演的结果和经过Tikh正则化的反演结果;图5A、图5B分别显示了中间变量B经过反演的结果和经过Tikh正则化的反演结果;图6A、图6B分别显示了各向异性梯度Bani经过反演的结果和经过Tikh正则化的反演结果;从该图中的可以看出,利用本发明的反演方法进行反演计算后,能较好的反演出地下各向异性的情况,说明本发明的反演方法有效性显著。
在介绍了本发明示例性实施方式的方法之后,接下来,参考图7对本发明示例性实施方式的基于稀疏表征的地震各向异性梯度反演系统进行介绍。
基于稀疏表征的地震各向异性梯度反演系统的实施可以参见上述方法的实施,重复之处不再赘述。以下所使用的术语“模块”或者“单元”,可以是实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现,但是硬件,或者软件和硬件的组合的实现也是可能并被构想的。
基于同一发明构思,本发明还提出了一种基于稀疏表征的地震各向异性梯度反演系统,如图7所示,该系统包括:
数据采集模块701,用于采集测井数据;
参数计算模块702,用于根据所述测井数据得到中间变量及各向异性梯度;
字典训练模块703,用于根据所述中间变量建立训练集,通过K-SVD方法从所述训练集中学习中间变量对应的稀疏字典;
插值处理模块704,用于对所述中间变量进行插值处理,建立各向异性介质中各向异性梯度反演的初始模型;
目标函数建立模块705,用于根据各向异性梯度理论和稀疏表征理论建立目标函数,其中,假设所述目标函数在同一工区中且地层在横向上具有连续性,将所述中间变量对应的稀疏字典上的稀疏表达作为约束条件;
反演计算模块706,用于将所述初始模型作为初始值,利用所述目标函数进行反演计算,得到中间变量的反演结果,并根据所述各向异性梯度与中间变量的关系得到各向异性梯度的反演结果。
在一实施例中,所述参数计算模块702具体用于:
根据所述测井数据得到中间变量A、B、C、D及各向异性梯度Bani,其中,
Figure BDA0002609256360000141
Figure BDA0002609256360000142
C=Bani cos 2φiso
D=Bani sin 2φiso
Figure BDA0002609256360000143
Figure BDA0002609256360000144
其中,A、B、C、D为中间变量;Bani为各向异性梯度;Biso为各向同性梯度;φiso为入射角;ΔZ为上下界面的纵波阻抗之差;
Figure BDA0002609256360000151
为上下界面的纵波阻抗的均值;Δδ、Δγ分别为上下界面的各项异性参数之差;
Figure BDA0002609256360000152
为上下界面的纵波速度平均值;
Figure BDA0002609256360000153
为上下界面的横波速度平均值;Δα为上下界面的纵波速度之差;ΔG为上下界面的横波剪切模量之差;
Figure BDA0002609256360000154
为上下界面的横波剪切模量平均值。
在一实施例中,所述字典训练模块703具体用于:
根据所述中间变量A、B、C、D建立训练集,其中训练集中包含多个训练样本;
通过K-SVD方法从所述训练集中学习中间变量对应的稀疏字典,供反演计算使用。
在一实施例中,反演计算模块706具体用于:
将初始模型中的中间变量进行分块处理,得到各参数的N个长度为M的训练样本集;
设置目标函数中的正则化参数以及最大迭代次数;
以所述中间变量对应的稀疏字典作为中间变量对应的初始字典,根据所述初始字典确定对应的初始化的稀疏系数矩阵;
根据正则化参数、最大迭代次数及目标函数,基于正则化方法反演得到中间变量的结果;
将稀疏字典中的原子固定,采用OMP方法对训练样本集进行稀疏编码,更新初始化的稀疏系数矩阵;
根据更新后的稀疏系数矩阵及中间变量的结果,重构出对应反演层段的变量值,得到中间变量的反演结果,并根据所述各向异性梯度与中间变量的关系得到各向异性梯度的反演结果。
应当注意,尽管在上文详细描述中提及了基于稀疏表征的地震各向异性梯度反演系统的若干模块,但是这种划分仅仅是示例性的并非强制性的。实际上,根据本发明的实施方式,上文描述的两个或更多模块的特征和功能可以在一个模块中具体化。反之,上文描述的一个模块的特征和功能可以进一步划分为由多个模块来具体化。
基于前述发明构思,如图8所示,本发明还提出了一种计算机设备800,包括存储器810、处理器820及存储在存储器810上并可在处理器820上运行的计算机程序830,所述处理器820执行所述计算机程序830时实现前述基于稀疏表征的地震各向异性梯度反演方法。
基于前述发明构思,本发明提出了一种计算机可读存储介质,所述计算机可读存储介质存储有计算机程序,所述计算机程序被处理器执行时实现前述基于稀疏表征的地震各向异性梯度反演方法。
本发明提出的基于稀疏表征的地震各向异性梯度反演方法及系统,能够更准确的发现地下的各向异性情况,根据实际测井数据获得中间变量,基于K-SVD方法从实际数据中学习得中间变量参数各自的稀疏字典,并将反演的各参数在对应的稀疏字典上的稀疏表达作为约束条件,从而进行各参数的反演,所得到的反演结果更趋向于分辨率更高的测井数据,实现纵向分辨率的提高且反演结果更准确,为后续的储层预测提供有力的数据支持。
本领域内的技术人员应明白,本发明的实施例可提供为方法、系统、或计算机程序产品。因此,本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且,本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器,使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中,使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品,该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上,使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理,从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
最后应说明的是:以上所述实施例,仅为本发明的具体实施方式,用以说明本发明的技术方案,而非对其限制,本发明的保护范围并不局限于此,尽管参照前述实施例对本发明进行了详细的说明,本领域的普通技术人员应当理解:任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化,或者对其中部分技术特征进行等同替换;而这些修改、变化或者替换,并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围,都应涵盖在本发明的保护范围之内。因此,本发明的保护范围应所述以权利要求的保护范围为准。

Claims (10)

1.一种基于稀疏表征的地震各向异性梯度反演方法,其特征在于,该方法包括:
采集测井数据;
根据所述测井数据得到中间变量及各向异性梯度;
根据所述中间变量建立训练集,通过K-SVD方法从所述训练集中学习中间变量对应的稀疏字典;
对所述中间变量进行插值处理,建立各向异性介质中各向异性梯度反演的初始模型;
根据各向异性梯度理论和稀疏表征理论建立目标函数,其中,假设所述目标函数在同一工区中且地层在横向上具有连续性,将所述中间变量对应的稀疏字典上的稀疏表达作为约束条件;
将所述初始模型作为初始值,利用所述目标函数进行反演计算,得到中间变量的反演结果,并根据所述各向异性梯度与中间变量的关系得到各向异性梯度的反演结果。
2.根据权利要求1所述的基于稀疏表征的地震各向异性梯度反演方法,其特征在于,根据所述测井数据得到中间变量及各向异性梯度,包括:
根据所述测井数据得到中间变量A、B、C、D及各向异性梯度Bani,其中,
Figure FDA0002609256350000011
Figure FDA0002609256350000012
C=Banicos2φiso
D=Banisin2φiso
Figure FDA0002609256350000013
Figure FDA0002609256350000014
其中,A、B、C、D为中间变量;Bani为各向异性梯度;Biso为各向同性梯度;φiso为入射角;ΔZ为上下界面的纵波阻抗之差;
Figure FDA0002609256350000021
为上下界面的纵波阻抗的均值;Δδ、Δγ分别为上下界面的各项异性参数之差;
Figure FDA0002609256350000022
为上下界面的纵波速度平均值;
Figure FDA0002609256350000023
为上下界面的横波速度平均值;Δα为上下界面的纵波速度之差;ΔG为上下界面的横波剪切模量之差;
Figure FDA0002609256350000024
为上下界面的横波剪切模量平均值。
3.根据权利要求2所述的基于稀疏表征的地震各向异性梯度反演方法,其特征在于,根据所述中间变量建立训练集,通过K-SVD方法从所述训练集中学习中间变量对应的稀疏字典,包括:
根据所述中间变量A、B、C、D建立训练集,其中训练集中包含多个训练样本;
通过K-SVD方法从所述训练集中学习中间变量对应的稀疏字典,供反演计算使用。
4.根据权利要求3所述的基于稀疏表征的地震各向异性梯度反演方法,其特征在于,将所述初始模型作为初始值,利用所述目标函数进行反演计算,得到中间变量的反演结果,并根据所述各向异性梯度与中间变量的关系得到各向异性梯度的反演结果,包括:
将初始模型中的中间变量进行分块处理,得到各参数的N个长度为M的训练样本集;
设置目标函数中的正则化参数以及最大迭代次数;
以所述中间变量对应的稀疏字典作为中间变量对应的初始字典,根据所述初始字典确定对应的初始化的稀疏系数矩阵;
根据正则化参数、最大迭代次数及目标函数,基于正则化方法反演得到中间变量的结果;
将稀疏字典中的原子固定,采用OMP方法对训练样本集进行稀疏编码,更新初始化的稀疏系数矩阵;
根据更新后的稀疏系数矩阵及中间变量的结果,重构出对应反演层段的变量值,得到中间变量的反演结果,并根据所述各向异性梯度与中间变量的关系得到各向异性梯度的反演结果。
5.一种基于稀疏表征的地震各向异性梯度反演系统,其特征在于,该系统包括:
数据采集模块,用于采集测井数据;
参数计算模块,用于根据所述测井数据得到中间变量及各向异性梯度;
字典训练模块,用于根据所述中间变量建立训练集,通过K-SVD方法从所述训练集中学习中间变量对应的稀疏字典;
插值处理模块,用于对所述中间变量进行插值处理,建立各向异性介质中各向异性梯度反演的初始模型;
目标函数建立模块,用于根据各向异性梯度理论和稀疏表征理论建立目标函数,其中,假设所述目标函数在同一工区中且地层在横向上具有连续性,将所述中间变量对应的稀疏字典上的稀疏表达作为约束条件;
反演计算模块,用于将所述初始模型作为初始值,利用所述目标函数进行反演计算,得到中间变量的反演结果,并根据所述各向异性梯度与中间变量的关系得到各向异性梯度的反演结果。
6.根据权利要求5所述的基于稀疏表征的地震各向异性梯度反演系统,其特征在于,所述参数计算模块具体用于:
根据所述测井数据得到中间变量A、B、C、D及各向异性梯度Bani,其中,
Figure FDA0002609256350000031
Figure FDA0002609256350000032
C=Banicos2φiso
D=Banisin2φiso
Figure FDA0002609256350000033
Figure FDA0002609256350000034
其中,A、B、C、D为中间变量;Bani为各向异性梯度;Biso为各向同性梯度;φiso为入射角;ΔZ为上下界面的纵波阻抗之差;
Figure FDA0002609256350000035
为上下界面的纵波阻抗的均值;Δδ、Δγ分别为上下界面的各项异性参数之差;
Figure FDA0002609256350000036
为上下界面的纵波速度平均值;
Figure FDA0002609256350000037
为上下界面的横波速度平均值;Δα为上下界面的纵波速度之差;ΔG为上下界面的横波剪切模量之差;
Figure FDA0002609256350000038
为上下界面的横波剪切模量平均值。
7.根据权利要求6所述的基于稀疏表征的地震各向异性梯度反演系统,其特征在于,所述字典训练模块具体用于:
根据所述中间变量A、B、C、D建立训练集,其中训练集中包含多个训练样本;
通过K-SVD方法从所述训练集中学习中间变量对应的稀疏字典,供反演计算使用。
8.根据权利要求7所述的基于稀疏表征的地震各向异性梯度反演系统,其特征在于,反演计算模块具体用于:
将初始模型中的中间变量进行分块处理,得到各参数的N个长度为M的训练样本集;
设置目标函数中的正则化参数以及最大迭代次数;
以所述中间变量对应的稀疏字典作为中间变量对应的初始字典,根据所述初始字典确定对应的初始化的稀疏系数矩阵;
根据正则化参数、最大迭代次数及目标函数,基于正则化方法反演得到中间变量的结果;
将稀疏字典中的原子固定,采用OMP方法对训练样本集进行稀疏编码,更新初始化的稀疏系数矩阵;
根据更新后的稀疏系数矩阵及中间变量的结果,重构出对应反演层段的变量值,得到中间变量的反演结果,并根据所述各向异性梯度与中间变量的关系得到各向异性梯度的反演结果。
9.一种计算机设备,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,其特征在于,所述处理器执行所述计算机程序时实现权利要求1至4任一所述方法。
10.一种计算机可读存储介质,其特征在于,所述计算机可读存储介质存储有计算机程序,所述计算机程序被处理器执行时实现权利要求1至4任一所述方法。
CN202010748575.7A 2020-07-30 基于稀疏表征的地震各向异性梯度反演方法及系统 Active CN114063156B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010748575.7A CN114063156B (zh) 2020-07-30 基于稀疏表征的地震各向异性梯度反演方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010748575.7A CN114063156B (zh) 2020-07-30 基于稀疏表征的地震各向异性梯度反演方法及系统

Publications (2)

Publication Number Publication Date
CN114063156A true CN114063156A (zh) 2022-02-18
CN114063156B CN114063156B (zh) 2024-06-25

Family

ID=

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106815876A (zh) * 2016-12-30 2017-06-09 清华大学 图像稀疏表征多字典学习的联合优化训练方法
CN110516358A (zh) * 2019-08-28 2019-11-29 电子科技大学 一种基于稀疏表示的地震各向异性参数反演方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106815876A (zh) * 2016-12-30 2017-06-09 清华大学 图像稀疏表征多字典学习的联合优化训练方法
CN110516358A (zh) * 2019-08-28 2019-11-29 电子科技大学 一种基于稀疏表示的地震各向异性参数反演方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张真梁: "基于字典学习稀疏表征的地震各向异性反演方法研究", 《中国优秀硕士学位论文全文数据库 基础科学辑》, no. 7, pages 8 - 39 *
李春鹏 等: "裂缝型储层预测的各向异性梯度反演方法研究", 《石油物探》, vol. 56, no. 6, pages 835 - 840 *

Similar Documents

Publication Publication Date Title
US11693139B2 (en) Automated seismic interpretation-guided inversion
CA3122488C (en) Subsurface models with uncertainty quantification
CN106405651B (zh) 一种基于测井匹配的全波形反演初始速度模型构建方法
EP2020609A1 (en) Methods and apparatus for geophysical exploration via joint inversion
US11181653B2 (en) Reservoir characterization utilizing ReSampled seismic data
CN113341465B (zh) 方位各向异性介质的地应力预测方法、装置、介质及设备
CN110516358A (zh) 一种基于稀疏表示的地震各向异性参数反演方法
Mousavi et al. Applications of deep neural networks in exploration seismology: A technical survey
CN116047583A (zh) 基于深度卷积神经网络的自适应波阻抗反演方法及系统
CN115455828A (zh) 一种可解释的测井曲线补全方法
CN111025388A (zh) 一种多波联合的叠前波形反演方法
CN111948718B (zh) 页岩气储层总有机碳含量预测方法及装置
CN114063156B (zh) 基于稀疏表征的地震各向异性梯度反演方法及系统
CN114063156A (zh) 基于稀疏表征的地震各向异性梯度反演方法及系统
CN116088048A (zh) 包含不确定性分析的各向异性介质多参数全波形反演方法
CN113970789B (zh) 全波形反演方法、装置、存储介质及电子设备
Zhang Ensemble methods of data assimilation in porous media flow for non-Gaussian prior probability density
CN112363222A (zh) 叠后自适应宽带约束波阻抗反演方法及装置
Li et al. A characterization method for cavity Karst reservoir using local full-waveform inversion in frequency domain
Torres et al. A deep-learning inverse Hessian preconditioning for iterative least-squares migration
CN115480306B (zh) 一种基于深度学习的隧道全波形反演优化方法及系统
Kordjazi et al. The use of the spectral element method for modeling stress wave propagation in non-destructive testing applications for drilled shafts
Parasyris et al. Near surface full waveform inversion via deep learning for subsurface imaging
CN118191911A (zh) 基于半监督学习的地震数据处理方法、装置、设备及介质
CN115657131A (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