CN113589381B - 一种基于压缩感知的相位与反射系数同时反演方法 - Google Patents

一种基于压缩感知的相位与反射系数同时反演方法 Download PDF

Info

Publication number
CN113589381B
CN113589381B CN202110910292.2A CN202110910292A CN113589381B CN 113589381 B CN113589381 B CN 113589381B CN 202110910292 A CN202110910292 A CN 202110910292A CN 113589381 B CN113589381 B CN 113589381B
Authority
CN
China
Prior art keywords
phase
reflection coefficient
wavelet
matrix
formula
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
CN202110910292.2A
Other languages
English (en)
Other versions
CN113589381A (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.)
Exploration Branch China Petroleum & Chemical Co Rporation
Chengdu Univeristy of Technology
Original Assignee
Exploration Branch China Petroleum & Chemical Co Rporation
Chengdu Univeristy of Technology
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 Exploration Branch China Petroleum & Chemical Co Rporation, Chengdu Univeristy of Technology filed Critical Exploration Branch China Petroleum & Chemical Co Rporation
Priority to CN202110910292.2A priority Critical patent/CN113589381B/zh
Publication of CN113589381A publication Critical patent/CN113589381A/zh
Application granted granted Critical
Publication of CN113589381B publication Critical patent/CN113589381B/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/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • 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/624—Reservoir parameters
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

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

本发明公开了一种基于压缩感知的相位与反射系数同时反演方法，相位谱分解包括：地震记录s表示为：s＝Gm+n，其中n是数据噪声，通过建立第一范数L1和第二范数L2最小化约束条件可以求解m系数，min[||s‑Gm||2+λ||m||1]，通过求解m并将序列m排列为矩阵形式，即得到
Figure DDA0003203453140000011
剖面，此时得到了相位
Figure DDA0003203453140000012
再将
Figure DDA0003203453140000013
剖面按相位相加则得到反射系数R。本发明提高了地震数据相位获取的直观性和准确性能较为直观的呈现地震数据的相位，同时也能得到去除相位变化影响的反射系数序列，本发明既能获取不同时间采样点的真实相位值，还能通过向量求和得到反射系数。

Description

一种基于压缩感知的相位与反射系数同时反演方法
技术领域
本发明涉及地震勘测技术领域，特别涉及一种基于压缩感知的相位与反射系数同时反演方法。
背景技术
地震反射系数反演是在频率域实施的，可获得高分辨率的时间域反射系数信息，实现对小于调谐厚度的薄储层进行有效识别，以达到提高地震数据分辨率和提高储层预测精度的目的。
反射系数在地震勘探表征一种去掉地震波波形影响之后的结果，它用来更直观的显示地层的层数和位置，一般在地震勘探中为了与测井曲线结果更好地对应，会对反射系数进行积分得到地震波阻抗，用地震波阻抗来显示地层。通过反射系数和它得到的波阻抗，能反应地震勘探中的地下的地层分布情况。通过这个地层的分布情况来定量描述含油气储层的分布情况。反射系数的求取是一个很精细的科研工作，传统的相位获取方法都是频率域方法，主要以时频分析方法获取相位谱，由于时频方法获取频谱和相位谱都是与频率相关，所以很难利用时频相位谱变换为与时间相关的相位谱。现有反射系数的技术在地震勘探的现象中存在分辨率不高，连续性不好，以及出现虚假或者多余的反射系数，这些现象都会导致积分后得到的地震波阻抗不容易识别储层。
目前的技术仅是单纯地从地震资料中提取反射系数，一是提取的信息仅为相对反射系数，二是在提取反射系数的过程中，并没有考虑相位的在波形信号上的影响，容易出现反射系数的错误匹配。而这些原因之所以存在该缺陷，是因为在所求解的褶积模型没有考虑子波相位的变化问题，所以在使用优化算法求解后，结果并没有体现对相位的鲁棒性。
发明内容
本发明的目的在于解决上述现有技术存在的缺陷，从时间域匹配相位的思路入手，以变相位子波库和时域地震信号为出发点，提供一种基于压缩感知的相位与反射系数同时反演方法，该方法利用时域压缩感知反演算法，逐相位求解，呈现了正确的时间域相位分布，解决了以往相位获取方法的不直观性，不稳定性，实现了对实际工区地震信号的相位刻画，以解决上述背景技术中提出的问题。
为实现上述目的，本发明提供如下技术方案：
一种基于压缩感知的相位与反射系数同时反演方法，相位谱分解包括以下步骤：
步骤1：地震记录s(t)表示为式(1)中子波与反射系数褶积的形式：
s(t)＝w(t)*r(t) (1)
其中，w(t)代表子波库、r(t)代表反射系数，式(1)写为式(2)中的矩阵形式：
s＝Wr (2)
W代表子波库，s代表地震记录，r代表反射系数；
步骤2：利用BP基追踪算法求解式(2)中未知量r的方法是通过构造如下的误差约束公式进行求解：
Ξ＝||s-Wr||2+λ||r||1 (3)
式(3)中，Ξ为目标函数，下标1和2分别代表向量的L1范数和L2范数，λ为正则化参数因子；满足Ξ→min的r即为最终的解向量；
步骤3：构造相位子波库：
G＝[Wφ1(t) Wφ2(t) … … Wφn(t)] (4)
式(4)是包含不同相位子波矩阵，相位为φi，i＝1…n；每个Wφi都是一个子波对角矩阵；G代表混合相位的子波库，Wφ1(t)代表不同相位子波方阵。
进一步地，任何反射系数r(t)序列都写成下式：
Figure BDA0003203453120000031
rΦi(t)代表不同相位的反射系数向量。
进一步地，地震记录s(t)表示为(6)的形式：
Figure BDA0003203453120000032
相位子波库[Wφ1(t) Wφ2(t) … … Wφn(t)]用矩阵G表示，反射系数序列[rφ1(t)rφ2(t) … … rφn(t)]T用矩阵m表示，其中T表示序列的转置，则地震记录s表示为：
s＝Gm+n (7)
其中n是数据噪声。
进一步地，通过建立第一范数L1和第二范数L2最小化约束条件可以求解m系数，由(3)式写为：
min[||s-Gm||2+λ||m||1] (8)
通过求解m并将序列m排列为矩阵形式，即得到
Figure BDA0003203453120000033
剖面，此时得到了相位/>
Figure BDA0003203453120000034
再将
Figure BDA0003203453120000035
剖面按相位相加则得到反射系数R。
进一步地，根据
Figure BDA0003203453120000036
剖面将反射系数的振幅值映射为相位值，得到以相位值为强度的真实相位剖面。
与现有技术相比，本发明的有益效果是：
1、本发明提出的基于压缩感知的相位与反射系数同时反演方法，提高了地震数据相位获取的直观性和准确性，能较为直观的呈现地震数据的相位，并且由于在解卷积的过程中子波库按相位分布，得到的反射系数不包含相位影响，是更加纯粹的反射系数。
2、本发明提出的基于压缩感知的相位与反射系数同时反演方法，利用子波按相位排列得到按相位分布的反射系数，获取了无相位影响的反射系数和不同采样点的真实相位，称为一种基于压缩感知的相位与反射系数同时反演方法。该方法并未扩大子波库的规模，而只是以原有数据规模重排子波库，因而既保证了计算效率又在传统方法的基础上获得了更多的信息。
附图说明
图1为相位库示意图；
图2为本发明的地震数据相位分解示意图；
图3为本发明的地震数据相位与反射系数同时反演示意图；
图4为本发明的ZJ研究区域地震数据；
图5为本发明的基于稀疏反演全谱相位分解的
Figure BDA0003203453120000041
追踪反演反射系数；
图6为本发明的基于稀疏反演全谱相位分解的
Figure BDA0003203453120000042
追踪的反演道积分；
图7为本发明的基于稀疏反演全谱相位分解的
Figure BDA0003203453120000043
追踪的相位剖面。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
一种基于压缩感知的相位与反射系数同时反演方法，包括基追踪求取相位
Figure BDA0003203453120000044
和反射系数R。
Figure BDA0003203453120000045
追踪原理
步骤1：地震记录s(t)表示为式(1)中子波与反射系数褶积的形式：
s(t)＝w(t)*r(t) (1)
其中，w(t)代表子波库、r(t)代表反射系数，式(1)写为式(2)中的矩阵形式：
s＝Wr (2)
W代表子波库，s代表地震记录，r代表反射系数；
步骤2：利用BP基追踪算法求解式(2)中未知量r的方法是通过构造如下的误差约束公式进行求解：
Ξ＝||s-Wr||2+λ||r||1 (3)
式(3)中，Ξ为目标函数，下标1和2分别代表向量的L1范数和L2范数，λ为正则化参数因子；满足Ξ→min的r即为最终的解向量；
步骤3：构造相位子波库：
G＝[Wφ1(t) Wφ2(t) … … Wφn(t)] (4)
式(4)是包含不同相位子波矩阵，相位为φi，i＝1…n；每个Wφi都是一个子波对角矩阵；G代表混合相位的子波库，Wφ1(t)代表不同相位子波方阵。
任何反射系数r(t)序列都写成下式：
Figure BDA0003203453120000051
rΦi(t)代表不同相位的反射系数向量。
地震记录s(t)表示为(6)的形式：
Figure BDA0003203453120000052
相位子波库[Wφ1(t) Wφ2(t) … … Wφn(t)]用矩阵G表示，反射系数序列[rφ1(t)rφ2(t) … … rφn(t)]T用矩阵m表示，其中T表示序列的转置，则地震记录s表示为：
s＝Gm+n (7)
其中n是数据噪声。
通过建立第一范数L1和第二范数L2最小化约束条件可以求解m系数，由(3)式写为：
min[||s-Gm||2+λ||m||1] (8)
通过求解m并将序列m排列为矩阵形式，即得到
Figure BDA0003203453120000066
剖面，此时得到了相位/>
Figure BDA0003203453120000067
再将
Figure BDA0003203453120000068
剖面按相位相加则得到反射系数R，本方法先得到相位，再得到反射系数；而本发明由于是算法获取相位后，再对相位进行求和得到反射系数，从整体技术来看为反射系数和相位的同时反演。
模型试算
为了验证基于稀疏反演全谱分解的相位与反射系数同时反演的可行性，利用图2最左边事先设置好的模型，利用文中提出的方法对其进行求解，得到图2中间的
Figure BDA0003203453120000061
剖面，该剖面的横坐标为相位，纵坐标为时间，垂直于/>
Figure BDA0003203453120000062
平面的轴为振幅，进一步的将/>
Figure BDA0003203453120000063
平面每个点变换为相应相位的相位子波，则得到原模型记录的相位波分解平面图，如图2右图。图3为相位反射系数同时反演的示意图，最左边为模型，同样右文中提出方法得到/>
Figure BDA0003203453120000064
剖面，将该剖面按/>
Figure BDA0003203453120000065
轴相加，则得到反射系数，通过模型的求解和讨论，该方法实现了相位和反射系数的同时反演，证明了方法的有效性。
本方法实现了相位和反射系数的同时反演，解决了提取的信息仅为相对反射系数，求取结果单一以及在提取反射系数的过程中，并没有考虑相位的在波形信号上的影响，容易出现反射系数的错误匹配的问题。
反射系数作积分可以得到地震波阻抗，在物理意义上，相位与反射系数和地震波阻抗没什么关系，但是在求解方法和数值上，一列反射系数可以视作不同相位的反射系数的叠加，相当于对地震反射系数按相位进行了一种分解，比如也可以按不同的频率来分解反射系数，如果当不同频率反射系数加到一起后，那么这一列数据就具有了所有频率信息。所以相位在数值上与反射系数有叠加和分解的关系，反射系数在数值上与地震波阻抗有积分的关系。
由于目前的技术仅是单纯地从地震资料中提取反射系数，一是提取的信息仅为相对反射系数，二是在提取反射系数的过程中，并没有考虑相位的在波形信号上的影响，容易出现反射系数的错误匹配。而这些原因之所以存在该缺陷，是因为在所求解的褶积模型没有考虑子波相位的变化问题，所以在使用优化算法求解后，结果并没有体现对相位的鲁棒性，本发明基于此，为了解决提取信息单一，将常规的反射系数求解方法进行改进，使结果在常规方法只能求解反射系数的基础上，也得到了地震数据的相位分布；为了解决相位变化带来的波形变化产生的反射系数错误匹配，将子波库按相位来排列，得到按相位排列的反射系数，解决了相位变化产生的反演问题。
以下介绍应用实例
以某区ZJ实际数据为实例，对其进行稀疏反演全谱相位分解的
Figure BDA0003203453120000071
追踪反演的验证，图4为该区数据的实际地震剖面，图5为利用文中方法对图4剖面所求的反射系数，可以看到横向连续性和纵向分辨率都很好，基于该反射系数作道积分得到地震波阻抗剖面得到图6，可以看到该阻抗的分辨率高，连续性也较好，对于薄互层的识别有较好的效果，且与井数据吻合；通过对/>
Figure BDA0003203453120000072
剖面进行数值映射方法，将反射系数的振幅值映射为相位值，得到真实相位剖面如图7，该相位剖面的分辨率很高，横向连续性较好，且与井数据吻合良好，验证了该方法的有效性，印证了利用相位子波库，同时求取反射系数和相位的可行性，准确性。
综上所述：本发明基于压缩感知的相位与反射系数同时反演方法，从算法所用的子波库入手，引入更加完备的具有相位信息的子波库，从而得到按相位和时间排列分布的反射系数，即
Figure BDA0003203453120000073
剖面，同时得到了反射系数和相位信息，根据/>
Figure BDA0003203453120000081
剖面将反射系数的振幅值映射为相位值，得到以相位值为强度的真实相位剖面，在对地震数据求解分相位反射系数得到了时间与相位的反射系数剖面，能够精确刻画地震数据的相位分布，从而得到高分辨率的分相位反射系数。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (1)

1.一种基于压缩感知的相位与反射系数同时反演方法，其特征在于，相位谱分解包括以下步骤：
步骤1：地震记录s(t)表示为式(1)中子波与反射系数褶积的形式：
s(t)＝w(t)*r(t) (1)
其中，w(t)代表子波库、r(t)代表反射系数，式(1)写为式(2)中的矩阵形式：
s＝Wr (2)
W代表子波库，s代表地震记录，r代表反射系数；
步骤2：利用BP基追踪算法求解式(2)中未知量r的方法是通过构造如下的误差约束公式进行求解：
Ξ＝||s-Wr||2+λr||1 (3)
式(3)中，Ξ为目标函数，下标1和2分别代表向量的L1范数和L2范数，λ为正则化参数因子；满足Ξ→min的r即为最终的解向量；
步骤3：构造相位子波库：
G＝[Wφ1(t) Wφ2(t)……Wφn(t)] (4)
式(4)是包含不同相位子波矩阵，相位为φi，i＝1…n；每个Wφi都是一个子波对角矩阵；G代表混合相位的子波库，Wφ1(t)代表不同相位子波方阵；
任何反射系数r(t)序列都写成下式：
Figure FDA0004223842660000011
rΦi(t)代表不同相位的反射系数向量；
地震记录s(t)表示为(6)的形式：
Figure FDA0004223842660000021
相位子波库[Wφ1(t) Wφ2(t)……Wφn(t)]用矩阵G表示，反射系数序列[rφ1(t) rφ2(t)……rφn(t)]T用矩阵m表示，其中T表示序列的转置，则地震记录s表示为：
s＝Gm+n (7)
其中n是数据噪声；
通过建立第一范数L1和第二范数L2最小化约束条件可以求解m系数，由(3)式写为：
min[||s-Gm||2+λ||m||1] (8)
通过求解m并将序列m排列为矩阵形式，即得到
Figure FDA0004223842660000022
剖面，此时得到了相位/>
Figure FDA0004223842660000023
再将/>
Figure FDA0004223842660000024
剖面按相位相加则得到反射系数R；
根据
Figure FDA0004223842660000025
剖面将反射系数的振幅值映射为相位值，得到以相位值为强度的真实相位剖面。
CN202110910292.2A 2021-08-09 2021-08-09 一种基于压缩感知的相位与反射系数同时反演方法 Active CN113589381B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110910292.2A CN113589381B (zh) 2021-08-09 2021-08-09 一种基于压缩感知的相位与反射系数同时反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110910292.2A CN113589381B (zh) 2021-08-09 2021-08-09 一种基于压缩感知的相位与反射系数同时反演方法

Publications (2)

Publication Number Publication Date
CN113589381A CN113589381A (zh) 2021-11-02
CN113589381B true CN113589381B (zh) 2023-06-27

Family

ID=78256554

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110910292.2A Active CN113589381B (zh) 2021-08-09 2021-08-09 一种基于压缩感知的相位与反射系数同时反演方法

Country Status (1)

Country Link
CN (1) CN113589381B (zh)

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4608643B2 (ja) * 2003-01-17 2011-01-12 株式会社武田エンジニアリング・コンサルタント 地震の予知方法、地震の予知システム、地震の予知プログラム及び記録媒体
US9075159B2 (en) * 2011-06-08 2015-07-07 Chevron U.S.A., Inc. System and method for seismic data inversion
CN103376464B (zh) * 2012-04-13 2016-04-06 中国石油天然气集团公司 一种地层品质因子反演方法
CN103018775B (zh) * 2012-11-15 2016-05-11 中国石油天然气股份有限公司 基于相位分解的混合相位子波反褶积方法
CN103852788A (zh) * 2014-02-27 2014-06-11 中国海洋石油总公司 一种基于复地震道分解和重构的地震相位和频率校正方法
WO2016065356A1 (en) * 2014-10-24 2016-04-28 Ion Geophysical Corporation Methods for seismic inversion and related seismic data processing
CN108594304A (zh) * 2018-07-25 2018-09-28 中国石油化工股份有限公司胜利油田分公司勘探开发研究院 基于多极子库利用线性规划求解l1范数的阻抗反演方法
CN109946739A (zh) * 2019-03-15 2019-06-28 成都理工大学 一种基于压缩感知理论的地震剖面增强方法
CN112444867A (zh) * 2019-08-30 2021-03-05 中国石油化工股份有限公司 一种基于压缩感知的地震资料频带拓展处理方法
CN112526599B (zh) * 2019-09-17 2024-04-09 中国石油化工股份有限公司 基于加权l1范数稀疏准则的子波相位估计方法及系统
CN111208561B (zh) * 2020-01-07 2020-09-01 自然资源部第一海洋研究所 基于时变子波与曲波变换约束的地震声波阻抗反演方法
CN113031068B (zh) * 2021-02-24 2022-05-27 浙江大学 一种基于反射系数精确式的基追踪叠前地震反演方法

Also Published As

Publication number Publication date
CN113589381A (zh) 2021-11-02

Similar Documents

Publication Publication Date Title
CN107817527B (zh) 基于块稀疏压缩感知的沙漠地震勘探随机噪声压制方法
CN106291677B (zh) 一种基于匹配追踪方法的叠后声波阻抗反演方法
CN110471104B (zh) 基于智能特征学习的叠后地震反射模式识别方法
CN105353408B (zh) 一种基于匹配追踪的Wigner高阶谱地震信号谱分解方法
CN102928875B (zh) 基于分数阶傅里叶域的子波提取方法
CN105093315B (zh) 一种去除煤层强反射信号的方法
CN104730576A (zh) 基于Curvelet变换的地震信号去噪方法
CN106291682A (zh) 一种基于基追踪方法的叠后声波阻抗反演方法
CN107121701A (zh) 基于Shearlet变换的多分量地震数据Corssline方向波场重建方法
Zhao et al. Signal detection and enhancement for seismic crosscorrelation using the wavelet-domain Kalman filter
CN113589381B (zh) 一种基于压缩感知的相位与反射系数同时反演方法
CN112213782B (zh) 分相位地震数据的处理方法、装置和服务器
Lin et al. Time-frequency mixed domain multi-trace simultaneous inversion method
CN114624765B (zh) 一种相位域地震数据处理与重构方法、装置及可存储介质
Wu et al. The suppression of powerline noise for TEM with coded source based on independent component analysis
CN110045418B (zh) 一种点坝侧积体三维地震识别方法
CN115755172A (zh) 基于构造约束压缩感知的薄互层高分辨率处理与识别方法
CN114200522B (zh) 深度域地震子波提取方法和装置、存储介质及电子设备
CN102353991B (zh) 基于匹配地震子波的物理小波的地震瞬时频率分析方法
CN110673211B (zh) 一种基于测井与地震数据的品质因子建模方法
Steinmann et al. Machine learning analysis of seismograms reveals a continuous plumbing system evolution beneath the Klyuchevskoy volcano in Kamchatka, Russia
Henriques et al. Improving the analysis of well-logs by wavelet cross-correlation
CN106125138B (zh) 一种叠前crp道集拓频处理的方法
CN117348081A (zh) 一种地震信号高频扩展方法及装置
CN113325471B (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