CN116088049B - 基于子波变换的最小二乘逆时偏移地震成像方法及装置 - Google Patents

基于子波变换的最小二乘逆时偏移地震成像方法及装置 Download PDF

Info

Publication number
CN116088049B
CN116088049B CN202310363163.5A CN202310363163A CN116088049B CN 116088049 B CN116088049 B CN 116088049B CN 202310363163 A CN202310363163 A CN 202310363163A CN 116088049 B CN116088049 B CN 116088049B
Authority
CN
China
Prior art keywords
wavelet
gaussian
sound wave
gaussian wavelet
seismic imaging
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
CN202310363163.5A
Other languages
English (en)
Other versions
CN116088049A (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.)
Tsinghua University
Original Assignee
Tsinghua University
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 Tsinghua University filed Critical Tsinghua University
Priority to CN202310363163.5A priority Critical patent/CN116088049B/zh
Publication of CN116088049A publication Critical patent/CN116088049A/zh
Application granted granted Critical
Publication of CN116088049B publication Critical patent/CN116088049B/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/282—Application of seismic models, synthetic seismograms
    • 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/303—Analysis for determining velocity profiles or travel times
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • 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)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Remote Sensing (AREA)
  • General Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Theoretical Computer Science (AREA)
  • Biomedical Technology (AREA)
  • Molecular Biology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Health & Medical Sciences (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提出基于子波变换的最小二乘逆时偏移地震成像方法及装置，属于地震成像技术领域。其中，所述方法包括：通过向待探测地质发射探测声波，获取接收器处接收到的探测声波信号；将所述探测声波信号输入预设的空洞卷积网络，所述空洞卷积网络输出拟合得到的所述接收器处的高斯子波信号；利用所述高斯子波信号进行最小二乘逆时偏移地震成像，得到所述待探测地质的地震成像结果。本发明可提高现有最小二乘逆时偏移地震成像方法的成像精度和成像速度，具有很高的应用价值。

Description

基于子波变换的最小二乘逆时偏移地震成像方法及装置
技术领域
本发明属于地震成像技术领域，特别提出一种基于子波变换的最小二乘逆时偏移地震成像方法及装置。
背景技术
地震成像是地震勘探中的关键环节之一。地震成像的本质是一种通过检波仪接收的人工地震产生的反射波，推断地下介质的速度结构等性质的地球物理勘探技术。目前逆时偏移成像是效果最好的地震成像方法之一。若将基于波动方程和先验模型参数的地震传播过程线性化，则逆时偏移成像可以被看作该线性化正演算子的伴随算子。该伴随算子是一个从地震观测数据空间到速度模型空间的映射。在计算过程中，根据此伴随算子对应的成像条件求出地下反射层的速度结构。计算成像条件时，需要将翻转时间的接收信号当做波源，在先验速度模型下反向传播计算反传波场，“逆时”的名字也因此而来。不过由于该方法只做一次成像，因此成像质量较差，且只能表示反射面，无法生成准确的速度模型。
而最小二乘逆时偏移地震成像是一种改进的逆时偏移成像算法，相当于一种线性化的波形反演，旨在寻找在最小二乘意义上最能预测记录的地震数据的速度模型。在计算过程中，先对波方程的正演算子线性化，再将线性化算子作用在残差速度模型上拟合残差数据，求出的残差速度模型即为成像结果。相比于逆时偏移成像而言，最小二乘逆时偏移地震成像可以减少成像中的伪影，提高成像分辨率。不过目前由于最小二乘逆时偏移地震成像是一种迭代算法，所需的计算时间很长，成像速度慢，因此尚未在工业界被普遍应用。
发明内容
本发明的目的是为克服已有技术的不足之处，提出一种基于子波变换的最小二乘逆时偏移地震成像方法及装置。本发明具有成像精度高、成像速度快的特点，弥补了已有最小二乘逆时偏移地震成像的缺陷。
本发明第一方面实施例提出一种基于子波变换的最小二乘逆时偏移地震成像方法，包括：
通过向待探测地质发射探测声波，获取接收器处接收到的探测声波信号；
将所述探测声波信号输入预设的空洞卷积网络，所述空洞卷积网络输出拟合得到的所述接收器处的高斯子波信号；
利用所述高斯子波信号进行最小二乘逆时偏移地震成像，得到所述待探测地质的地震成像结果。
在本发明的一个具体实施例中，所述空洞卷积网络包括：依次连接的第一卷积层、多个空洞卷积层、第二卷积层和一个全连接层；
其中，所述第一卷积层包含卷积和参数化整流线性单元两种操作，所述第二卷积层包含卷积操作。
在本发明的一个具体实施例中，所述利用所述高斯子波信号进行最小二乘逆时偏移地震成像，包括：
地震成像的正问题表示为：
Figure SMS_1
其中，
Figure SMS_11
表示/>
Figure SMS_5
的残差地震数据，/>
Figure SMS_8
是/>
Figure SMS_14
的速度模型扰动，/>
Figure SMS_16
是/>
Figure SMS_17
的矩阵，/>
Figure SMS_18
是信号采样点数量，/>
Figure SMS_12
是速度模型的维度；/>
Figure SMS_13
是以先验速度模型/>
Figure SMS_6
为背景模型的所述空洞卷积网络输出的接收器处拟合的高斯子波信号/>
Figure SMS_10
与该接收器接收的高斯子波信号/>
Figure SMS_4
之差，即/>
Figure SMS_7
；/>
Figure SMS_9
是初始高斯子波波场/>
Figure SMS_15
在接收器处接收的信号，
Figure SMS_2
，/>
Figure SMS_3
由求解如下波方程得到：
Figure SMS_19
其中，
Figure SMS_22
表示在时刻/>
Figure SMS_23
处于位置/>
Figure SMS_26
的初始高斯子波波场，/>
Figure SMS_20
代表拉普拉斯算子，/>
Figure SMS_24
；/>
Figure SMS_27
分别表示位置/>
Figure SMS_30
在平行地表方向和垂直地表方向的坐标；/>
Figure SMS_21
表示在时刻/>
Figure SMS_25
处于位置/>
Figure SMS_28
的初始高斯子波波场；在时刻
Figure SMS_29
时，地质中任一位置的初始高斯子波波场为0；
Figure SMS_31
代表发射机的位置；/>
Figure SMS_32
是狄拉克函数，/>
Figure SMS_33
表示在位置/>
Figure SMS_34
处的点源，/>
Figure SMS_35
表示在时刻/>
Figure SMS_36
高斯波源发射的声波强度；/>
Figure SMS_37
代表探测总时长；
则逆时偏移地震成像问题为如下所示的最小二乘问题：
Figure SMS_38
求解所述最小二乘问题，得到
Figure SMS_39
的解即为最终的地震成像结果。
在本发明的一个具体实施例中，所述空洞卷积层共有11层，每个空洞卷积层的大小为
Figure SMS_40
，其中第i个空洞卷积层的扩张率为2i-1，i=1,2,…,11；
所述第一卷积层、所述第二卷积层和所述空洞卷积层的卷积核数量均为64个。
在本发明的一个具体实施例中，在所述将所述探测声波信号输入预设的空洞卷积网络之前，所述方法还包括：训练所述空洞卷积网络；
所述训练所述空洞卷积网络包括：
1）构建训练集，所述训练集中的训练样本由仿真中在同一接收器处收到的探测声波信号和高斯子波信号组成，所述探测声波信号和所述高斯子波信号分别通过预设的探测声波方程和高斯子波声波方程获取；
2）构建所述空洞卷积网络；
3）利用所述训练集训练所述空洞卷积网络，得到训练完毕的所述空洞卷积网络。
在本发明的一个具体实施例中，所述探测声波方程表达式如下：
Figure SMS_41
其中，
Figure SMS_60
表示声波在/>
Figure SMS_44
处的传播速度；/>
Figure SMS_55
表示在时刻/>
Figure SMS_53
处于位置/>
Figure SMS_57
的探测声波波场；/>
Figure SMS_59
代表拉普拉斯算子，/>
Figure SMS_61
，/>
Figure SMS_51
分别表示位置/>
Figure SMS_56
在平行地表方向和垂直地表方向的坐标；/>
Figure SMS_42
表示在时刻/>
Figure SMS_48
发射机发射的探测声波强度；/>
Figure SMS_50
代表发射机的位置；/>
Figure SMS_52
是狄拉克函数，/>
Figure SMS_54
用来表示在位置/>
Figure SMS_58
处的点源；/>
Figure SMS_45
代表探测总时长；/>
Figure SMS_49
表示在时刻/>
Figure SMS_46
处于位置/>
Figure SMS_47
的探测声波波场；在时刻/>
Figure SMS_43
时，地质中任一位置探测声波的波场为0；
所述高斯子波声波方程表达式如下：
Figure SMS_62
其中，
Figure SMS_65
表示在时刻/>
Figure SMS_68
处于位置/>
Figure SMS_71
的高斯子波波场；
Figure SMS_63
，/>
Figure SMS_67
表示在时刻/>
Figure SMS_69
高斯波源发射的声波强度；/>
Figure SMS_72
表示在时刻/>
Figure SMS_64
处于位置/>
Figure SMS_66
的高斯子波波场；在时刻/>
Figure SMS_70
时，地质中任一位置高斯子波的波场为0。
在本发明的一个具体实施例中，所述高斯子波的频率大于等于40Hz。
本发明第二方面实施例提出一种基于子波变换的最小二乘逆时偏移地震成像装置，包括：
探测声波获取模块，用于通过向待探测地质发射探测声波，获取接收器处接收到的探测声波信号；
高斯子波拟合模块，用于将所述探测声波信号输入预设的空洞卷积网络，所述空洞卷积网络输出拟合得到的所述接收器处的高斯子波信号；
成像模块，用于利用所述高斯子波信号进行最小二乘逆时偏移地震成像，得到所述待探测地质的地震成像结果。
本发明第三方面实施例提出一种电子设备，包括：
至少一个处理器；以及，与所述至少一个处理器通信连接的存储器；
其中，所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被设置为用于执行上述一种基于子波变换的最小二乘逆时偏移地震成像方法。
本发明第四方面实施例提出一种计算机可读存储介质，所述计算机可读存储介质存储计算机指令，所述计算机指令用于使所述计算机执行上述一种基于子波变换的最小二乘逆时偏移地震成像方法。
本发明的特点及有益效果：
本发明设计了一种用高斯子波信号加速最小二乘逆时偏移地震成像的技术，具有成像速度快、成像精度高的特点。本发明将特定波形的声波发射到地质中后用接收器接收地震探测声波信号，然后将此探测声波信号用预设的神经网络转换成高斯子波信号，使用高斯子波信号会大幅加速最小二乘偏逆时偏移成像。本发明与传统的最小二乘逆时偏移地震成像方法相比，成像精度高且成像速度可提高30%左右，大大节约了成像的资源及时间。
本发明的应用场景广泛，在各种地质结构的速度成像问题上均可使用。
附图说明
本发明上述的和/或附加的方面和优点从下面结合附图对实施例的描述中将变得明显和容易理解，其中：
图1是本发明实施例中一种基于子波变换的最小二乘逆时偏移地震成像方法的整体流程图；
图2是本发明一个具体实施例中用于将探测声波信号转换为高斯子波信号的空洞卷积网络的结构示意图；
图3是本发明一个具体实施例中空洞卷积核的示意图。
具体实施方式
本发明提出一种基于子波变换的最小二乘逆时偏移地震成像方法及装置，下面结合附图和具体实施例进一步详细说明如下。
本发明第一方面实施例提出一种基于子波变换的最小二乘逆时偏移地震成像方法，包括：
通过向待探测地质发射探测声波，获取接收器处接收到的探测声波信号；
将所述探测声波信号输入预设的空洞卷积网络，所述空洞卷积网络输出拟合得到的所述接收器处的高斯子波信号；
利用所述高斯子波信号进行最小二乘逆时偏移地震成像，得到所述待探测地质的地震成像结果。
在本发明的一个具体实施例中，所述基于子波变换的最小二乘逆时偏移地震成像方法，该方法整体流程如图1所示，包括训练阶段和测试阶段，具体步骤如下：
1）训练阶段。
1-1）获取训练集。
1-1-1）分别构建用于描述探测声波和高斯子波在地质中的传播过程的声波方程。
本实施例中，单个训练样本由任一时刻任一位置的探测声波信号和对应的高斯子波信号组成。
具体地，探测声波在地质中的传播由只含速度参数的声波方程决定。在本发明一个具体实施例中，探测声波在地质中的传播过程可由如下声波方程决定：
Figure SMS_73
其中，
Figure SMS_84
表示声波在/>
Figure SMS_76
处的传播速度，由预设的速度模型确定，在本发明一个具体实施例中设置为Marmousi速度模型；/>
Figure SMS_79
表示在时刻/>
Figure SMS_89
处于位置/>
Figure SMS_91
的探测声波波场；
Figure SMS_88
代表拉普拉斯算子，/>
Figure SMS_92
，其中/>
Figure SMS_85
分别表示位置/>
Figure SMS_90
在平行地表方向和垂直地表方向的坐标；/>
Figure SMS_77
表示在时刻/>
Figure SMS_80
发射机发射的探测声波强度，在本发明一个具体实施例中设置为频率为20 Hz的瑞克子波；/>
Figure SMS_74
代表发射机的位置；/>
Figure SMS_78
是狄拉克函数，/>
Figure SMS_82
表示在位置/>
Figure SMS_86
处的点源。/>
Figure SMS_81
代表探测总时长。/>
Figure SMS_83
表示在时刻/>
Figure SMS_87
处于位置/>
Figure SMS_93
的探测声波波场。在时刻/>
Figure SMS_75
时，地质中任一位置探测声波的波场为0，即所探测地质中没有声波存在。
本实施例中，由于已知发射机位置、探测声波波形和声波速度，故求解上述带有初值条件的声波方程式（1），即得到在任一时刻t任一位置
Figure SMS_94
处的探测声波波场/>
Figure SMS_95
，从而可得对应/>
Figure SMS_96
的接收器处收到的探测声波信号/>
Figure SMS_97
，其中/>
Figure SMS_98
代表接收器的位置，/>
Figure SMS_99
即为训练集中的输入数据。在本发明一个具体实施例中探测声波波场设置为瑞克子波波场，则接收器处收到的为瑞克子波信号。
高斯子波在地质中的传播过程可由如下声波方程决定：
Figure SMS_100
(2)
其中，
Figure SMS_103
表示在时刻/>
Figure SMS_104
处于位置/>
Figure SMS_108
的高斯子波波场；/>
Figure SMS_101
代表拉普拉斯算子，
Figure SMS_105
，/>
Figure SMS_109
表示在时刻/>
Figure SMS_110
高斯波源发射的声波强度，本实施例中高斯子波的频率设置为40Hz及以上，在本发明一个具体实施例中设置为频率为50 Hz的高斯子波。/>
Figure SMS_102
表示在时刻/>
Figure SMS_106
处于位置/>
Figure SMS_107
的高斯子波波场。在时刻/>
Figure SMS_111
时，地质中任一位置高斯子波的波场为0。
接收器处收到时刻t位置
Figure SMS_112
处的的高斯子波信号记为/>
Figure SMS_113
。/>
Figure SMS_114
即为训练集中的输出标签。相同位置的接收器接收到的/>
Figure SMS_115
和/>
Figure SMS_116
组成一个训练样本。
1-1-2）根据步骤1-1-1）的声波方程，通过构建仿真地震场景获取训练样本。
本实施例中，首先随机选取5-10个速度模型作为训练速度模型，在本发明一个具体实施例中，选取Marmousi模型的5个切片作为训练速度模型，由于本实施例仿真的是海上地震数据，因此在每个训练速度模型的顶部添加了一个深度为190m的水层，陆地数据则不需要添加水层。选取不同速度模型的原因是增强神经网络的泛化性。其次，本发明实施例中采用神经网络的输入的震源波形应是探测声波的波形，输出的震源波形应是高频的高斯子波。在本发明一个具体实施例中，输入的震源波形应是以20Hz为主频的高通的瑞克小波，输出的源震源波形是以50Hz为中心的高斯小波；训练时的地震采集观测系统与测试时的地震采集观测系统相同。
在本发明一个具体实施例中，在仿真中，将231个震源均匀地分布在地表深20m处；共设置920个接收器，空间采样间隔为10m，所有接收器均匀地分布在地表深30m处；数据采集时，在每个训练速度模型下，令每个震源单独起震一次，然后通过步骤1-1-1）建立的声波方程构分别获取探测声波波场数据和对应的高斯子波波场数据以构建训练样本，本实施例中，每次采集波场数据包含2048个时间采样点，时间采样间隔为2.16ms；将每次地震位置相邻3-5个接收器接收到的探测声波信号和对应高斯子波信号组成一个训练样本。在本发明一个具体实施例中，把5个相邻位置的接收器接收的探测声波信号数据和对应高斯子波信号数据作为一个训练样本，因此训练样本中训练数据的大小为
Figure SMS_117
，训练样本的个数为/>
Figure SMS_118
。
本实施例中，按照上述设定的训练速度模型、地震场景和数据采集形式，采用有限差分法求解步骤1-1-1）的声波方程，生成二维合成地震数据作为训练样本，所有训练样本组成训练集。
1-2）构建用于将探测声波信号转换为高斯子波信号的空洞卷积网络。
图2为本发明一个具体实施例中用于将探测声波信号转换为高斯子波信号的空洞卷积网络的结构示意图。如图2所示，该网络由依次连接的第一卷积层、多个空洞卷积层、第二卷积层和一个全连接层组成。
具体地，所述空洞卷积网络的第一层为第一卷积层，所述第一卷积层的输入是探测声波信号，该卷积层包含两种类型的操作：卷积和参数化整流线性单元(ReLU)。其中，第一卷积层的卷积滤波器可以写成：
Figure SMS_119
其中，
Figure SMS_120
是该卷积层的输入在/>
Figure SMS_121
位置的取值，/>
Figure SMS_122
是卷积算子，/>
Figure SMS_123
表示滤波器权重，/>
Figure SMS_124
是卷积核大小。
参数化整流线性单元可以表述为：
Figure SMS_125
其中，
Figure SMS_126
是所述空洞卷积网络在ReLU之前的输入。第一卷积层的主要作用是提取不同层次的特征。在本发明的一个具体实施例中，第一卷积层的输入的大小为/>
Figure SMS_127
的训练数据，输出为/>
Figure SMS_128
的向量。本实施例中，所述第一卷积层的输出进入依次连接的一系列的/>
Figure SMS_129
的空洞卷积层。所述空洞卷积层的空洞卷积核是在常规的卷积核元素之间插入零，从而跳过输入状态的一些点。本实施例中，空洞卷积的表达式为：
Figure SMS_130
其中，
Figure SMS_131
是下取整算子，/>
Figure SMS_132
是设定的扩张率。在本发明一个具体实施例中，输入数据在时间方向上的长度为2048，因此本发明的一个具体实施例使用了11个依次连接的空洞卷积层，其扩张率分别为/>
Figure SMS_133
，11个依次连接的空洞卷积层的感受野大小为2048，可以聚合全时范围的信息。图3显示了本发明一个具体实施例中在时间维度（图3中横向）上扩张率为/>
Figure SMS_134
的空洞卷积核，/>
Figure SMS_135
表示卷积核的元素，扩张率控制了卷积核中相邻元素间的距离。对于空间维度（图3中纵向），空洞卷积核采取常规卷积的形式。空洞卷积层的主要作用是扩大神经网络的感受野。在本发明的一个具体实施例中，空洞卷积层的输出为/>
Figure SMS_136
的向量。最后一个空洞卷积层的输出结果输入至第二卷积层，所述第二卷积层只有卷积操作而没有激活函数。在本发明的一个具体实施例中，第二卷积层的输出为/>
Figure SMS_137
的向量。本实施例中，第一卷积层、第二卷积层和空洞卷积层均包含64个卷积核。第二卷积层的输出结果输入至全连接层，全连接层数学上可以表述为：
Figure SMS_138
其中，z表示当前输入，
Figure SMS_139
和/>
Figure SMS_140
表示待优化参数。
令
Figure SMS_141
和/>
Figure SMS_142
表示输入的探测声波信号的第/>
Figure SMS_143
个输入数据和输出的高斯子波信号的第/>
Figure SMS_144
个标签，则损失函数/>
Figure SMS_145
可以表述为：
Figure SMS_146
其中，
Figure SMS_147
表示训练样本的总数，/>
Figure SMS_148
表示逐点的/>
Figure SMS_149
损失函数。
1-3）利用步骤1-1）得到的训练集训练步骤1-2）建立的空洞卷积网络，得到训练完毕的空洞卷积网络；
在本发明一个具体实施例中，以5个相邻接收器的信号组成一个训练样本，因此神经网络的输入和输出都是
Figure SMS_150
。网络训练时使用Adam优化器更新网络权重；网络初始化采用正交初始化方法；学习率为循环余弦退火学习率，本实施例的初始学习率为/>
Figure SMS_151
，并在32代中衰减到/>
Figure SMS_152
，共训练80代；训练时的批大小（即每次调整参数前所选取的样本数量）为32。
训练完成后，得到训练完毕的空洞卷积网络。
2）测试阶段。
2-1）利用超声波发射机向待探测地质发射探测声波，获取接收器处接收到的探测声波信号。
本实施例中，所述超声波发射机所发射的探测声波具有给定的波形和频率，其中所述超声波发射机的型号无特殊要求；探测声波的波形和频率由具体的发射机给出，无特殊要求。
总体而言，测试阶段的实验场景及发射器波形与训练阶段一致。具体来说，在本发明的一个具体实施例中，探测声波为频率20 Hz的瑞克子波信号。探测声波发射机的间距一般为40m-100m，均匀布置在地表或者海平面下10m-30m，在本发明一个具体实施例中发射机间距设置为40m，布置于待探测的海平面下20m，在探测过程中在不同的时间点发射声波，以避免互相干扰。发射器个数没有限制，在本发明的一个具体实施例中，使用一个发射机在不同的时间点发射探测声波。接收器的间距一般为5m-30m，均匀布置在地表或者海平面下10m-30m，接收器个数根据探测区域大小以及接收器间距而定，在本发明一个具体实施例中接收器的间距设置为10m，布置于待探测的海平面下30m，个数为921个。接收器收到的信号间隔一般认为可以非常小，在本发明一个具体实施例中收到信号的间隔设置为2.16ms。
2-2）将步骤2-1）获取的探测声波信号输入步骤1）训练完毕的空洞卷积网络，所述网络输出对应输入的探测声波信号的拟合得到的在对应探测器处的高斯子波信号，表达式如下：
Figure SMS_153
其中，
Figure SMS_154
为拟合得到的高斯子波信号。
2-3）根据步骤2-2）得到的高斯子波信号，根据高斯子波对应的声波方程进行最小二乘逆时偏移地震成像。
本实施例中，地震成像的正问题可以表示为：
Figure SMS_155
,
其中，
Figure SMS_167
表示/>
Figure SMS_159
的残差地震数据，/>
Figure SMS_162
是/>
Figure SMS_170
的速度模型扰动，/>
Figure SMS_172
是/>
Figure SMS_171
的矩阵，该矩阵是取决于数据采集形式、震源波形和已知速度模型的线性化波动方程正算子，这里/>
Figure SMS_174
是信号采样点数量，/>
Figure SMS_165
是速度模型的维度，在本发明的一个具体实例中分别为
Figure SMS_168
，/>
Figure SMS_156
。/>
Figure SMS_161
是以先验速度模型/>
Figure SMS_157
为背景模型的接收器处拟合的高斯子波信号/>
Figure SMS_160
与该接收器接收的高斯子波信号/>
Figure SMS_164
之差，即/>
Figure SMS_169
。其中/>
Figure SMS_158
是初始高斯子波波场/>
Figure SMS_163
在接收器处接收的信号，/>
Figure SMS_166
，/>
Figure SMS_173
由求解如下波方程给出：
Figure SMS_175
其中，
Figure SMS_177
表示在时刻/>
Figure SMS_179
处于位置/>
Figure SMS_182
的初始高斯子波波场，
Figure SMS_178
。/>
Figure SMS_180
表示在时刻/>
Figure SMS_181
处于位置/>
Figure SMS_183
的初始高斯子波波场。在时刻/>
Figure SMS_176
时，地质中任一位置的初始高斯子波波场为0。
则逆时偏移地震成像问题可以写成如下最小二乘问题：
Figure SMS_184
,
该问题可使用已有的高效优化算法求解，得到的结果
Figure SMS_185
即为最终的地震成像结果。在本发明一个具体实施例中使用共轭梯度法求解。
为实现上述实施例，本发明第二方面实施例提出一种基于子波变换的最小二乘逆时偏移地震成像装置，包括：
探测声波获取模块，用于通过向待探测地质发射探测声波，获取接收器处接收到的探测声波信号；
高斯子波拟合模块，用于将所述探测声波信号输入预设的空洞卷积网络，所述空洞卷积网络输出拟合得到的所述接收器处的高斯子波信号；
成像模块，用于利用所述高斯子波信号进行最小二乘逆时偏移地震成像，得到所述待探测地质的地震成像结果。
需要说明的是，前述对一种基于子波变换的最小二乘逆时偏移地震成像方法的实施例解释说明也适用于本实施例的一种基于子波变换的最小二乘逆时偏移地震成像装置，在此不再赘述。根据本发明实施例提出的一种基于子波变换的最小二乘逆时偏移地震成像装置，通过向待探测地质发射探测声波，获取接收器处接收到的探测声波信号；将所述探测声波信号输入预设的空洞卷积网络，所述空洞卷积网络输出拟合得到的所述接收器处的高斯子波信号；利用所述高斯子波信号进行最小二乘逆时偏移地震成像，得到所述待探测地质的地震成像结果。由此可实现提高现有最小二乘逆时偏移地震成像方法的成像精度和成像速度，具有很高的应用价值。
为实现上述实施例，本发明第三方面实施例提出一种电子设备，包括：
至少一个处理器；以及，与所述至少一个处理器通信连接的存储器；
其中，所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被设置为用于执行上述一种基于子波变换的最小二乘逆时偏移地震成像方法。
为实现上述实施例，本发明第四方面实施例提出一种计算机可读存储介质，所述计算机可读存储介质存储计算机指令，所述计算机指令用于使所述计算机执行上述基于子波变换的最小二乘逆时偏移地震成像方法。
需要说明的是，本公开上述的计算机可读介质可以是计算机可读信号介质或者计算机可读存储介质或者是上述两者的任意组合。计算机可读存储介质例如可以是——但不限于——电、磁、光、电磁、红外线、或半导体的系统、装置或器件，或者任意以上的组合。计算机可读存储介质的更具体的例子可以包括但不限于：具有一个或多个导线的电连接、便携式计算机磁盘、硬盘、随机访问存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、光纤、便携式紧凑磁盘只读存储器(CD-ROM)、光存储器件、磁存储器件、或者上述的任意合适的组合。在本公开中，计算机可读存储介质可以是任何包含或存储程序的有形介质，该程序可以被指令执行系统、装置或者器件使用或者与其结合使用。而在本公开中，计算机可读信号介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了计算机可读的程序代码。这种传播的数据信号可以采用多种形式，包括但不限于电磁信号、光信号或上述的任意合适的组合。计算机可读信号介质还可以是计算机可读存储介质以外的任何计算机可读介质，该计算机可读信号介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。计算机可读介质上包含的程序代码可以用任何适当的介质传输，包括但不限于：电线、光缆、RF(射频)等等，或者上述的任意合适的组合。
上述计算机可读介质可以是上述电子设备中所包含的；也可以是单独存在，而未装配入该电子设备中。上述计算机可读介质承载有一个或者多个程序，当上述一个或者多个程序被该电子设备执行时，使得该电子设备执行上述实施例的一种基于子波变换的最小二乘逆时偏移地震成像方法。
可以以一种或多种程序设计语言或其组合来编写用于执行本公开的操作的计算机程序代码，上述程序设计语言包括面向对象的程序设计语言—诸如Java、Smalltalk、C++，还包括常规的过程式程序设计语言—诸如“C”语言或类似的程序设计语言。程序代码可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络——包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本申请的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本申请的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。
流程图中或在此以其他方式描述的任何过程或方法描述可以被理解为，表示包括一个或更多个用于实现特定逻辑功能或过程的步骤的可执行指令的代码的模块、片段或部分，并且本申请的优选实施方式的范围包括另外的实现，其中可以不按所示出或讨论的顺序，包括根据所涉及的功能按基本同时的方式或按相反的顺序，来执行功能，这应被本申请的实施例所属技术领域的技术人员所理解。
在流程图中表示或在此以其他方式描述的逻辑和/或步骤，例如，可以被认为是用于实现逻辑功能的可执行指令的定序列表，可以具体实现在任何计算机可读介质中，以供指令执行系统、装置或设备(如基于计算机的系统、包括处理器的系统或其他可以从指令执行系统、装置或设备取指令并执行指令的系统)使用，或结合这些指令执行系统、装置或设备而使用。就本说明书而言，"计算机可读介质"可以是任何可以包含、存储、通信、传播或传输程序以供指令执行系统、装置或设备或结合这些指令执行系统、装置或设备而使用的装置。计算机可读介质的更具体的示例(非穷尽性列表)包括以下：具有一个或多个布线的电连接部(电子装置)，便携式计算机盘盒(磁装置)，随机存取存储器(RAM)，只读存储器(ROM)，可擦除可编辑只读存储器(EPROM或闪速存储器)，光纤装置，以及便携式光盘只读存储器(CDROM)。另外，计算机可读介质甚至可以是可在其上打印程序的纸或其他合适的介质，因为可以例如通过对纸或其他介质进行光学扫描，接着进行编辑、解译或必要时以其他合适方式进行处理来以电子方式获得程序，然后将其存储在计算机存储器中。
应当理解，本申请的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行系统执行的软件或固件来实现。例如，如果用硬件来实现，和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列(PGA)，现场可编程门阵列(FPGA)等。
本技术领域的普通技术人员可以理解实现上述实施例方法携带的全部或部分步骤是可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，该程序在执行时，包括方法实施例的步骤之一或其组合。
此外，在本申请各个实施例中的各功能单元可以集成在一个处理模块中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个模块中。上述集成的模块既可以采用硬件的形式实现，也可以采用软件功能模块的形式实现。集成的模块如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。
上述提到的存储介质可以是只读存储器，磁盘或光盘等。尽管上面已经示出和描述了本申请的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本申请的限制，本领域的普通技术人员在本申请的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (9)

1.一种基于子波变换的最小二乘逆时偏移地震成像方法，其特征在于，包括：
通过向待探测地质发射探测声波，获取接收器处接收到的探测声波信号；
将所述探测声波信号输入预设的空洞卷积网络，所述空洞卷积网络输出拟合得到的所述接收器处的高斯子波信号；
利用所述高斯子波信号进行最小二乘逆时偏移地震成像，得到所述待探测地质的地震成像结果；
其中，所述利用所述高斯子波信号进行最小二乘逆时偏移地震成像，包括：
地震成像的正问题表示为：
Figure QLYQS_1
其中，δd表示T×1的残差地震数据，δm是M×1的速度模型扰动，
Figure QLYQS_3
是T×M的矩阵，T是信号采样点数量，M是速度模型的维度；δd是以先验速度模型c0(x)为背景模型的所述空洞卷积网络输出的接收器处拟合的高斯子波信号/>
Figure QLYQS_6
与该接收器接收的高斯子波信号/>
Figure QLYQS_9
之差，即/>
Figure QLYQS_4
Figure QLYQS_5
是初始高斯子波波场/>
Figure QLYQS_7
在接收器处接收的信号，/>
Figure QLYQS_8
Figure QLYQS_2
由求解如下波方程得到：
Figure QLYQS_10
Figure QLYQS_11
其中，
Figure QLYQS_12
表示在时刻t处于位置x的初始高斯子波波场，Δ代表拉普拉斯算子，
Figure QLYQS_13
x1，x2分别表示位置x在平行地表方向和垂直地表方向的坐标；/>
Figure QLYQS_14
表示在时刻0处于位置x的初始高斯子波波场；在时刻0时，地质中任一位置的初始高斯子波波场为0；
xs代表发射机的位置；δ是狄拉克函数，δ(x-xs)表示在位置xs处的点源，fG(t)表示在时刻t高斯波源发射的声波强度；T代表探测总时长；
则逆时偏移地震成像问题为如下所示的最小二乘问题：
Figure QLYQS_15
求解所述最小二乘问题，得到δm的解即为最终的地震成像结果。
2.根据权利要求1所述的方法，其特征在于，所述空洞卷积网络包括：依次连接的第一卷积层、多个空洞卷积层、第二卷积层和一个全连接层；
其中，所述第一卷积层包含卷积和参数化整流线性单元两种操作，所述第二卷积层包含卷积操作。
3.根据权利要求2所述的方法，其特征在于，所述空洞卷积层共有11层，每个空洞卷积层的大小为2×3，其中第i个空洞卷积层的扩张率为2i-1，i＝1，2，…，11；
所述第一卷积层、所述第二卷积层和所述空洞卷积层的卷积核数量均为64个。
4.根据权利要求2所述的方法，其特征在于，在所述将所述探测声波信号输入预设的空洞卷积网络之前，所述方法还包括：训练所述空洞卷积网络；
所述训练所述空洞卷积网络包括：
1)构建训练集，所述训练集中的训练样本由仿真中在同一接收器处收到的探测声波信号和高斯子波信号组成，所述探测声波信号和所述高斯子波信号分别通过预设的探测声波方程和高斯子波声波方程获取；
2)构建所述空洞卷积网络；
3)利用所述训练集训练所述空洞卷积网络，得到训练完毕的所述空洞卷积网络。
5.根据权利要求4所述的方法，其特征在于，所述探测声波方程表达式如下：
Figure QLYQS_16
Figure QLYQS_17
其中，c(x)表示声波在x处的传播速度；uR(x，t)表示在时刻t处于位置x的探测声波波场；Δ代表拉普拉斯算子，
Figure QLYQS_18
x＝(x1，x2)，x1，x2分别表示位置x在平行地表方向和垂直地表方向的坐标；fR(t)表示在时刻t发射机发射的探测声波强度；xs代表发射机的位置；δ是狄拉克函数，δ(x-xs)用来表示在位置xs处的点源；T代表探测总时长；uR(x，0)表示在时刻0处于位置x的探测声波波场；在时刻0时，地质中任一位置探测声波的波场为0；
所述高斯子波声波方程表达式如下：
Figure QLYQS_19
Figure QLYQS_20
其中，uG(x，t)表示在时刻t处于位置x的高斯子波波场；
Figure QLYQS_21
fG(t)表示在时刻t高斯波源发射的声波强度；uG(x，0)表示在时刻0处于位置x的高斯子波波场；在时刻0时，地质中任一位置高斯子波的波场为0。
6.根据权利要求4所述的方法，其特征在于，所述高斯子波的频率大于等于40Hz。
7.一种基于子波变换的最小二乘逆时偏移地震成像装置，其特征在于，包括：
探测声波获取模块，用于通过向待探测地质发射探测声波，获取接收器处接收到的探测声波信号；
高斯子波拟合模块，用于将所述探测声波信号输入预设的空洞卷积网络，所述空洞卷积网络输出拟合得到的所述接收器处的高斯子波信号；
成像模块，用于利用所述高斯子波信号进行最小二乘逆时偏移地震成像，得到所述待探测地质的地震成像结果；
其中，所述利用所述高斯子波信号进行最小二乘逆时偏移地震成像，包括：
地震成像的正问题表示为：
Figure QLYQS_22
其中，δd表示T×1的残差地震数据，δm是M×1的速度模型扰动，
Figure QLYQS_25
是T×M的矩阵，T是信号采样点数量，M是速度模型的维度；δd是以先验速度模型c0(x)为背景模型的所述空洞卷积网络输出的接收器处拟合的高斯子波信号/>
Figure QLYQS_26
与该接收器接收的高斯子波信号/>
Figure QLYQS_29
之差，即/>
Figure QLYQS_24
Figure QLYQS_27
是初始高斯子波波场/>
Figure QLYQS_28
在接收器处接收的信号，/>
Figure QLYQS_30
Figure QLYQS_23
由求解如下波方程得到：
Figure QLYQS_31
Figure QLYQS_32
其中，
Figure QLYQS_33
表示在时刻t处于位置x的初始高斯子波波场，Δ代表拉普拉斯算子，
Figure QLYQS_34
x＝(x1，x2)；x1，x2分别表示位置x在平行地表方向和垂直地表方向的坐标；/>
Figure QLYQS_35
表示在时刻0处于位置x的初始高斯子波波场；在时刻0时，地质中任一位置的初始高斯子波波场为0；
xs代表发射机的位置；δ是狄拉克函数，δ(x-xs)表示在位置xs处的点源，fG(t)表示在时刻t高斯波源发射的声波强度；T代表探测总时长；
则逆时偏移地震成像问题为如下所示的最小二乘问题：
Figure QLYQS_36
求解所述最小二乘问题，得到δm的解即为最终的地震成像结果。
8.一种电子设备，其特征在于，包括：
至少一个处理器；以及，与所述至少一个处理器通信连接的存储器；
其中，所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被设置为用于执行上述权利要求1-6任一项所述的方法。
9.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储计算机指令，所述计算机指令用于使所述计算机执行权利要求1-6任一项所述的方法。
CN202310363163.5A 2023-04-07 2023-04-07 基于子波变换的最小二乘逆时偏移地震成像方法及装置 Active CN116088049B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310363163.5A CN116088049B (zh) 2023-04-07 2023-04-07 基于子波变换的最小二乘逆时偏移地震成像方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310363163.5A CN116088049B (zh) 2023-04-07 2023-04-07 基于子波变换的最小二乘逆时偏移地震成像方法及装置

Publications (2)

Publication Number Publication Date
CN116088049A CN116088049A (zh) 2023-05-09
CN116088049B true CN116088049B (zh) 2023-06-20

Family

ID=86204802

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310363163.5A Active CN116088049B (zh) 2023-04-07 2023-04-07 基于子波变换的最小二乘逆时偏移地震成像方法及装置

Country Status (1)

Country Link
CN (1) CN116088049B (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101334482A (zh) * 2008-08-05 2008-12-31 中国海洋石油总公司 一种预测地震波中的多次波和一次波信号的方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2019046550A1 (en) * 2017-09-01 2019-03-07 The Trustees Of Princeton University QUANTITATIVE ULTRASOUND IMAGING BASED ON SEISMIC COMPLETE WAVEFORM REVERSAL
CN110857998B (zh) * 2018-08-23 2021-11-05 中国石油化工股份有限公司 一种基于lowrank有限差分的弹性逆时偏移方法及系统
CN111290019B (zh) * 2020-03-16 2021-04-13 中国海洋大学 一种应用于最小二乘逆时偏移的l-bfgs初始矩阵求取方法
CN112083482B (zh) * 2020-08-06 2021-11-19 西安交通大学 基于模型驱动深度学习的地震超分辨反演方法
US11604299B2 (en) * 2020-12-18 2023-03-14 Pgs Geophysical As Mixed-phase source wavelet estimation from recorded seismic data
CN113688696B (zh) * 2021-08-04 2023-07-18 南京信息工程大学 一种超高分遥感影像震害建筑物检测方法
CN114428324B (zh) * 2022-04-06 2022-06-28 中国石油大学(华东) 叠前高角度快速傅里叶变换地震成像方法、系统、设备
CN115079264A (zh) * 2022-06-28 2022-09-20 东北石油大学 一种不依赖子波的粘声最小二乘逆时偏移成像方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101334482A (zh) * 2008-08-05 2008-12-31 中国海洋石油总公司 一种预测地震波中的多次波和一次波信号的方法

Also Published As

Publication number Publication date
CN116088049A (zh) 2023-05-09

Similar Documents

Publication Publication Date Title
US10802167B2 (en) Seismic acquisition method and apparatus
KR101092668B1 (ko) 파형 역산을 이용한 지하 구조의 영상화 장치와 방법
US8737165B2 (en) Interferometric seismic data processing for a towed marine survey
US9091787B2 (en) Separation of simultaneous source data
US11029432B2 (en) De-aliased source separation method
US10215869B2 (en) System and method of estimating anisotropy properties of geological formations using a self-adjoint pseudoacoustic wave propagator
KR20110057124A (ko) 지진 표면파들의 파형들을 사용하는 토양 특성들의 추정
US10877175B2 (en) Seismic acquisition geometry full-waveform inversion
US11467307B2 (en) Methods and data processing apparatus for deblending seismic data
CN103630931A (zh) 从近场测量和建模假想特征计算假想源特征的方法和系统
Symes Mathematics of reflection seismology
CN103901473B (zh) 一种基于非高斯性最大化的双检信号上下行波场分离方法
CN111665556B (zh) 地层声波传播速度模型构建方法
Prior et al. Characterization of the acoustic output of single marine-seismic airguns and clusters: The Svein Vaage dataset
Aulanier et al. Time-angle sensitivity kernels for sound-speed perturbations in a shallow ocean
CN116088049B (zh) 基于子波变换的最小二乘逆时偏移地震成像方法及装置
Wei et al. A new P-wave reconstruction method for VSP data using conditional generative adversarial network
CN115170428A (zh) 一种声波远探测成像图的降噪方法
AU2014374370A1 (en) System and method of mitigating instabilities in a pseudoacoustic wave propagator
EP3090276A1 (en) System and method of mitigating instabilities in a pseudoacoustic wave propagator
US9784869B2 (en) Noise models by selection of transform coefficients
Zhou et al. Flattening the seismic data for optimal noise attenuation
US11686872B2 (en) Attenuation of guided waves using polarization filtering
CN112114361B (zh) 一种地下浅层复杂空间中震动场时空层析成像方法
CN111694052B (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