基于物理约束神经网络的地震速度反演方法

基于物理约束神经网络的地震速度反演方法
1.1、一种基于物理约束神经网络的地震速度反演方法，其特征在于：地震速度反演方法包括以下步骤：
2.步骤一、训练数据集准备
3.速度反演任务建立由地震数据到速度模型的非线性映射建模，需要大量成对的地震数据-速度模型样本，但真实速度模型往往很难获得；采用人工生成二维速度模型(简单层状模型、断层模型和盐丘模型等)的方法扩充数据得到速度模型数据集d_vel；
4.人工生成速度模型过程中，通过正弦、余弦函数以及线性函数的叠加生成层数在layer_m～layer_n之间，且层速度在velocity_m～velocity_n范围内的简单层状速度模型，速度模型尺寸为velocity_len
×
velocity_wid，横向和纵向的网格大小均为grid_size；在此基础上生成断层长度、位置随机及倾斜角度在angle_m～angle_n的断层模型；同时将大小为salt_dome_size个网格点、形状随机、速度为salt_dome_vel的盐丘随机加入断层模型中得到盐丘体模型；
5.获得足够的速度模型后，采用有限差分正演的方法生成对应的理论地震数据集d_sei；对于震源和观测系统设置，采用地表放炮和地表观测方式，检波器在地表水平间隔每detector_interval_size放置一个(即地表每个网格点均有)，采样时间间隔为t，总采集时间为t，即接收到的地震道集记录有t/t个数据；震源子波采用主频为frequency_size的雷克子波，共shot_size炮，震源水平间隔为source_interval_size，地表逐点激发；
6.为了模拟真实地震数据含有噪声和缺失道的影响，在理论地震数据集d_sei的基础上随机对dataset_percent_a的地震数据加入强度不等的高斯噪声，噪声标准差系数σ最大为noise_dev，另对dataset_percent_b的地震数据进行随机缺失地震道处理，最大随机缺失比为dele_percent，这样得到观测地震数据d_sei_obs，扩充数据集的同时可以提高网络模型的鲁棒性；
7.步骤二：先验知识处理和引入
8.对于震源位置和检波器位置信息采用one-hot编码的方式引入到网络中；由于震源和检波器位置信息是地震数据和速度模型关联的特征，设置两行编码向量，一行表示震源位置信息，一行表示检波器位置信息，将其拼接到地震数据集d_sei_obs中，记为数据集d_sei_obs_pos；然后设置大小和速度模型一致的两个二维矩阵，通过one-hot编码方式一个反映震源位置，一个反映检波器位置,得到数据集d_vel_pos；测井数据作为地震勘探中可靠的先验知识，可以作为反演方法的地质先验信息约束；通过well_num口竖井的测井资料获取到速度模型数据集d_vel中第col1列到第col
well_num
列，并使用最近邻插值算法和高斯模糊操作生成粗糙的地下初始速度模型作为初始速度先验知识,记为数据集d_init；
9.步骤三：网络输入与标签设定
10.以步骤二中得到的加入震源和检波器位置信息的地震数据d_sei_obs_pos作为第一个编码器的输入，反映震源和检波器位置的数据集d_vel_pos和生成粗糙的地下初始速度模型数据集d_init作为第二个编码器的输入，以对应的速度模型数据集d_vel作为标签，为训练神经网络模型g提供数据；
11.步骤四：基于物理先验约束的速度反演网络g构建
12.网络模型g分为基于数据驱动的速度反演子网和基于rnn正演约束子网组成；基于数据驱动的速度反演子网由一个双编码器和一个解码器组成；基于rnn正演约束子网结合全波形反演理论和rnn网络，实现了物理规律与神经网络的有效结合；
13.网络模型设定为
14.(1)、数据驱动的速度反演子网gd15.速度反演子网gd的两个编码器分别输入加入震源和检波器位置信息的地震数据和生成粗糙的地下初始速度模型；通过测井先验提供大尺度波速模型，而地震数据的波形信息通过另一组编码器网络提供小尺度特征信息；通过两个编码网络将地震数据特征压缩至一维向量形式，格式为[bs，c，1，1]，其中bs是批处理大小、c为特征通道数，然后将提取的特征输入解码器得到网络的反演结果；
[0016]
(2)、rnn正演约束子网g
p
[0017]
将基于数据驱动的速度反演子网得到的反演结果作为rnn正演约束子网的输入，结合全波形反演理论和rnn网络处理时序数据的特征，将以波动方程驱动的波场延拓规律与rnn网络结合，通过卷积方式计算波场拉普拉斯值并利用隐藏状态层保存上个时刻的波场信息u，实现地震正演的过程，并结合神经网络自动微分算法和adam最优化算法实现自动求取梯度并更新参数模型；
[0018]
步骤五：损失函数设定
[0019]
整体网络的损失函数为：
[0020]
l＝l(x，y)+l(d
true
，d
pre
)
[0021]
x为速度模型数据集d_vel，y为速度反演子网的输出，d
true
和d
pre
分别为理论地震数据d_sei和rnn正演约束子网得到的预测地震数据；
[0022]
其中l(x，y)由均方误差函数和结构相似性函数(structural similarity，ssim)共同衡量x和y的差异，具体形式为：
[0023][0024][0025]
其中，μ
x
、μy表示x和y内样本元素的平均值，表示其方差，σ
xy
表示其协方差，c1和c2表示计算稳定性常数，n表示批样本数量；ssim是测量两个样本之间相似性的指标，主要关注边缘和纹理相似性，ssim通常取值范围为[0，1]，为1时表示两个样本完全相似，因此其取值范围从0开始越接近1越好，为了满足神经网络最小化损失函数的优化目标，采用负数形式满足要求；
[0026]
l(d
true
，d
pre
)采用mse损失函数衡量模拟地震数据和观测地震数据的误差，结果为：
[0027][0028]
步骤六：预训练数据驱动的速度反演子网
[0029]
考虑到rnn网络内存消耗较大，且训练两个串联网络调参难度较大，基于迁移学习的思想，使用上述步骤二得到的加入先验知识的数据集预训练速度反演子网，将预训练好的网络模型保存；
[0030]
步骤七：结合rnn正演约束网络联合训练
[0031]
将预训练的速度反演子网和rnn正演约束子网联合起来训练迭代，实现基于物理规律和先验知识约束的地震速度反演网络；该网络将反演网络输出的预测速度模型作为rnn正演约束网络的可训练参数输入使两个网络联通起来，在训练中经过rnn正演约束网络得到的预测地震数据和理论地震数据计算损失，实现对速度反演子网的梯度计算和参数更新，训练结束后，以反演子网的输出预测结果作为整体网络的最终速度反演结果。
[0032]

[0033]
本发明属于地球科学技术与人工智能交叉领域，具体涉及一种基于物理约束神经网络的地震速度反演方法。


背景技术：

[0034]
石油产业是关系到民生和国家安全的重要产业，石油与天然气是国家发展的重要战略资源，我国石油进口量逐年增加，国内原油的供需缺口也在逐年加大，中国己经成为了世界上最大的能源进口国和消费国。寻找地下油气资源，提高国内的油气勘探开发能力具有重要的战略意义。地震勘探技术是解决油气勘探问题效果最好的一种方法。油气地震勘探的根本任务是根据观测到的各种信息研究和提取有关地下介质的物性参数，如速度、密度等，并对储层的含油气性做出评价。完成这一任务有正演和反演两种途径，它们也是弹性动力学研究的两个基本方面。正演是在给定震源和介质特性时研究地震波的传播规律，而反演则是根据各种地球物理观测数据推测地球内部的结构、形态及物质成分，定量计算各种相关的地球物理参数，按照模型介质的复杂程度，可将物理参数分为速度、各向异性和吸收衰减参数，对这些参数的定量分析可以大大提高我们对储层的认识。地震波速度参数贯穿于地震数据采集、处理和解释的整个过程，是地震勘探中最重要的参数之一，对速度参数的反演是地震反演的主要工作之一。
[0035]
目前，利用接收到的地震数据对地下进行速度反演的方法主要有速度分析、层析成像和全波形反演三类。其中，速度分析和层析成像方法只利用到地震波场的运动学特征(走时)，这从理论上限制了这两类方法进行速度反演的精度；而全波形反演不仅利用了波场的运动学特征，还利用了波场的动力学特征，是目前地震勘探领域公认的速度反演精度最高的方法，但仍然面临许多问题。限制全波形反演技术广泛应用的问题主要有地震数据低频缺失问题、初始模型严重依赖性和计算效率问题等。
[0036]
深度学习是机器学习和人工智能研究的最新趋势之一。它也是当今最流行的科学研究趋势之一。深度学习方法为计算机视觉和机器学习带来了革命性的进步。新的深度学习技术正在不断诞生，超越最先进的机器学习甚至是现有的深度学习技术。近年来，全世界在这一领域取得了许多重大突破。深度学习是一种特征学习方法，把原始数据通过非线性的模型转变成为更高层次的、更加抽象的表达，通过足够多的转换的组合学习复杂的映射
关系。将深度学习技术引入到地震油气勘探速度反演领域之中,利用数据驱动的方法来对多模态数据进行深入挖掘,以提取大数据中隐含的可学习的复杂结构和规律，在一定程度上解决传统速度参数反演的瓶颈问题，进一步提高业界对高精度地震波成像和油藏描述的预期。
[0037]
深度学习为速度建模提供了一种高效且可以准确利用大量地震数据的手段。2018年araya-polo等提出了一种基于深度神经网络(dnn)的二维速度建模方法，其将事先提取的特征信息作为网络输入，速度模型作为网络输出，通过神经网络近似两者间的非线性映射。2019年zhang等开发一种基于实时数据驱动技术的速度建模对抗网络(velocitygan)，直接从原始地震波形数据生成高质量的速度图像，且结果优于物理驱动的波形反演方法，提高了速度模型精度。2020年yuan等在已知初始速度和数据的情况下，用全卷积神经网络(fcn)实现了时移数据与目标速度变化之间的映射。2021年孟涛等研究了一种基于深度学习优化的全波形反演方法，在反演过程中利用深度学习中的adam优化算法迭代更新速度模型参数，得到了较好的反演结果。但是，地震勘探速度反演领域的应用场景相对于传统深度学习领域的应用场景更具复杂性和多解性，使用纯数据驱动的深度学习方法与预期目标仍有差距。深度学习速度建模方法通过神经网络框架对数据空间进行全局搜索，克服了线性反演方法中初始模型依赖的限制，相比于传统方法计算效率更高，是一种非常有应用前景的地震速度反演方法。但是纯数据驱动的深度学习方法受数据质量、数量影响较大，在反演精度和泛化上仍有较大提升空间。
[0038]
上述可知，将基于物理模型驱动的速度反演方法与基于数据驱动的深度学习方法结合起来，构成基于数据与模型双驱动的速度反演方法。基于深度学习方法解决速度反演中的问题，结合数据与模型双驱动的思想，通过地震勘探方法提供更多的特征工程信息，来降低深度神经网络模型学习难度以及提高泛化能力和结果的精度。


技术实现要素：

[0039]
本发明的目的是提供一种结合数据与模型双驱动的神经网络地震速度反演方法与流程，结合测井数据先验和注意力机制提高网络反演效果，然后基于循环神经网络实现对声波波动方程求解的地震正演网络，利用并行计算的优势实现正演过程的快速计算，并利用深度学习中的自动微分算法完成梯度计算，实现对反演网络的梯度传导完成物理约束条件，最后结合两个网络结构构建出最终基于波场延拓规律约束的地震速度反演网络，提高地震速度反演网络模型的精度与泛化能力。
[0040]
本发明采用的技术方案为：基于物理约束神经网络的地震速度反演方法，地震速度反演方法包括以下步骤：
[0041]
步骤一、训练数据集准备
[0042]
速度反演任务建立由地震数据到速度模型的非线性映射建模，需要大量成对的地震数据-速度模型样本，但真实速度模型往往很难获得；采用人工生成二维速度模型(简单层状模型、断层模型和盐丘模型等)的方法扩充数据得到速度模型数据集d_vel；
[0043]
人工生成速度模型过程中，通过正弦、余弦函数以及线性函数的叠加生成层数在layer_m～layer_n之间，且层速度在velocity_m～velocity_n范围内的简单层状速度模型，速度模型尺寸为velocity_len
×
velocity_wid，横向和纵向的网格大小均为grid_
size；在此基础上生成断层长度、位置随机及倾斜角度在angle_m～angle_n的断层模型；同时将大小为salt_dome_size个网格点、形状随机、速度为salt_dome_vel的盐丘随机加入断层模型中得到盐丘体模型；
[0044]
获得足够的速度模型后，采用有限差分正演的方法生成对应的理论地震数据集d_sei；对于震源和观测系统设置，采用地表放炮和地表观测方式，检波器在地表水平间隔每detector_interval_size放置一个(即地表每个网格点均有)，采样时间间隔为t，总采集时间为t，即接收到的地震道集记录有t/t个数据；震源子波采用主频为frequency_size的雷克子波，共shot_size炮，震源水平间隔为source_interval_size，地表逐点激发；
[0045]
为了模拟真实地震数据含有噪声和缺失道的影响，在理论地震数据集d_sei的基础上随机对dataset_percent_a的地震数据加入强度不等的高斯噪声，噪声标准差系数σ最大为noise_dev，另对dataset_percent_b的地震数据进行随机缺失地震道处理，最大随机缺失比为dele_percent，这样得到观测地震数据d_sei_obs，扩充数据集的同时可以提高网络模型的鲁棒性；
[0046]
步骤二：先验知识处理和引入
[0047]
对于震源位置和检波器位置信息采用one-hot编码的方式引入到网络中；由于震源和检波器位置信息是地震数据和速度模型关联的特征，设置两行编码向量，一行表示震源位置信息，一行表示检波器位置信息，将其拼接到地震数据集d_sei_obs中，记为数据集d_sei_obs_pos；然后设置大小和速度模型一致的两个二维矩阵，通过one-hot编码方式一个反映震源位置，一个反映检波器位置,得到数据集d_vel_pos；测井数据作为地震勘探中可靠的先验知识，可以作为反演方法的地质先验信息约束；通过well_num口竖井的测井资料获取到速度模型数据集d_vel中第col1列到第col
well_num
列，并使用最近邻插值算法和高斯模糊操作生成粗糙的地下初始速度模型作为初始速度先验知识,记为数据集d_init；
[0048]
步骤三：网络输入与标签设定
[0049]
以步骤二中得到的加入震源和检波器位置信息的地震数据d_sei_obs_pos作为第一个编码器的输入，反映震源和检波器位置的数据集d_vel_pos和生成粗糙的地下初始速度模型数据集d_init作为第二个编码器的输入，以对应的速度模型数据集d_vel作为标签，为训练神经网络模型g提供数据；
[0050]
步骤四：基于物理先验约束的速度反演网络g构建
[0051]
网络模型g分为基于数据驱动的速度反演子网和基于rnn正演约束子网组成；基于数据驱动的速度反演子网由一个双编码器和一个解码器组成；基于rnn正演约束子网结合全波形反演理论和rnn网络，实现了物理规律与神经网络的有效结合；
[0052]
网络模型设定为
[0053]
(1)、数据驱动的速度反演子网gd[0054]
速度反演子网gd的两个编码器分别输入加入震源和检波器位置信息的地震数据和生成粗糙的地下初始速度模型；通过测井先验提供大尺度波速模型，而地震数据的波形信息通过另一组编码器网络提供小尺度特征信息；通过两个编码网络将地震数据特征压缩至一维向量形式，格式为[bs，c，1，1]，其中bs是批处理大小、c为特征通道数，然后将提取的特征输入解码器得到网络的反演结果；
[0055]
(2)、rnn正演约束子网g
p
[0056]
将基于数据驱动的速度反演子网得到的反演结果作为rnn正演约束子网的输入，结合全波形反演理论和rnn网络处理时序数据的特征，将以波动方程驱动的波场延拓规律与rnn网络结合，通过卷积方式计算波场拉普拉斯值并利用隐藏状态层保存上个时刻的波场信息u，实现地震正演的过程，并结合神经网络自动微分算法和adam最优化算法实现自动求取梯度并更新参数模型；
[0057]
步骤五：损失函数设定
[0058]
整体网络的损失函数为：
[0059]
l＝l(x，y)+l(d
true
，d
pre
)
[0060]
x为速度模型数据集d_vel，y为速度反演子网的输出，d
true
和d
pre
分别为理论地震数据d_sei和rnn正演约束子网得到的预测地震数据；
[0061]
其中l(x，y)由均方误差函数和结构相似性函数(structural similarity，ssim)共同衡量x和y的差异，具体形式为：
[0062][0063][0064]
其中，μ
x
、μy表示x和y内样本元素的平均值，表示其方差，σ
xy
表示其协方差，c1和c2表示计算稳定性常数，n表示批样本数量；ssim是测量两个样本之间相似性的指标，主要关注边缘和纹理相似性，ssim通常取值范围为[0，1]，为1时表示两个样本完全相似，因此其取值范围从0开始越接近1越好，为了满足神经网络最小化损失函数的优化目标，采用负数形式满足要求；
[0065]
l(d
true
，d
pre
)采用mse损失函数衡量模拟地震数据和观测地震数据的误差，结果为：
[0066][0067]
步骤六：预训练数据驱动的速度反演子网
[0068]
考虑到rnn网络内存消耗较大，且训练两个串联网络调参难度较大，基于迁移学习的思想，使用上述步骤二得到的加入先验知识的数据集预训练速度反演子网，将预训练好的网络模型保存；
[0069]
步骤七：结合rnn正演约束网络联合训练
[0070]
将预训练的速度反演子网和rnn正演约束子网联合起来训练迭代，实现基于物理规律和先验知识约束的地震速度反演网络；该网络将反演网络输出的预测速度模型作为rnn正演约束网络的可训练参数输入使两个网络联通起来，在训练中经过rnn正演约束网络得到的预测地震数据和理论地震数据计算损失，实现对速度反演子网的梯度计算和参数更新，训练结束后，以反演子网的输出预测结果作为整体网络的最终速度反演结果。
[0071]
本发明的有益效果：提供了一种结合数据与模型双驱动的神经网络地震速度反演方法与流程，为了降低反演问题的多解性引入测井、观测系统和震源位置的先验知识加快模型的收敛，然后结合注意力机制的u-net网络加强对先验知识的利用；同时又考虑到地震
数据与速度模型的弱空间对应关系和地震数据中反射波信息丰富但能量较弱的特点通过深层卷积将特征大小提取到1
×
1并结合全连接层利用到尽可能多的特征进行空间映射，实现了基于数据驱动的速度反演网络并证明了其有效性；结合全波形反演理论和rnn网络处理时序数据的特征，将以波动方程驱动的波场延拓规律与rnn网络结合，通过卷积方式计算波场拉普拉斯值并利用隐藏状态层保存上个时刻的波场信息u,实现地震正演的过程，提高了计算效率；并结合神经网络自动微分算法和adam最优化算法实现自动求取梯度并更新参数模型，代替了传统全波形反演方法，实现了物理规律与神经网络的有效结合；通过将数据驱动的神经网络和物理规律驱动的rnn正演网络结合起来实现基于数据和模型双驱动的速度反演网络，并结合迁移学习的思想先预训练数据驱动的速度反演网络，然后结合rnn正演约束网络进行再训练，最终网络在较少样本的情况下实现高精度的速度反演。与现有同类技术相比本发明在结果精度和模型泛化性上均有一定的优势。
附图说明
[0072]
图1是实施例一的流程图；
[0073]
图2是实施例一的基于物理先验约束的速度反演网络模型图；
[0074]
图3是实施例二中训练指标ssim曲线图；
[0075]
图4是实施例二中训练指标mse曲线图；
[0076]
图5是实施例二中不同缺失道断层测试结果对比图；
[0077]
图6是实施例二中不同缺失道盐丘测试结果对比图；
[0078]
图7是实施例二中不同噪声等级断层测试结果对比图；
[0079]
图8是实施例二中不同噪声等级盐丘测试结果对比图；
[0080]
图9是实施例二中同类算法对比结果图；
[0081]
图10是实施例二中同类算法在断层模型第70列速度值对比图。
具体实施方式
[0082]
实施例一
[0083]
参照图1和图2，一种基于物理约束神经网络的地震速度反演方法，地震速度反演方法包括以下步骤：
[0084]
步骤一、训练数据集准备
[0085]
速度反演任务建立由地震数据到速度模型的非线性映射建模，需要大量成对的地震数据-速度模型样本，但真实速度模型往往很难获得；采用人工生成二维速度模型(简单层状模型、断层模型和盐丘模型等)的方法扩充数据得到速度模型数据集d_vel；
[0086]
人工生成速度模型过程中，首先通过正弦、余弦函数以及各种线性函数的叠加生成层数在4层到7层之间且层速度在1500～4500m/s范围内的简单层状速度模型，速度模型尺寸为1000m
×
640m，横向和纵向的网格大小均为10m；在此基础上生成断层长度、位置随机且倾斜角度在30
°
～150
°
的断层模型；同时将大小为20
×
20个网格点、形状随机、速度为5km/s的盐丘随机加入上述速度模型中得到盐丘体模型；
[0087]
获得足够的速度模型后，采用有限差分正演的方法生成对应的理论地震数据集d_sei；首先对于震源和观测系统设置，本发明采用地表放炮和地表观测方式，检波器在地表
水平间隔每10m放置一个(即地表每个网格点均有)，采样时间间隔为0.001s，总采集时间为1.024s，即接收到的地震道集记录有1024个数据；震源子波采用主频为20hz的雷克子波，共15炮，震源水平间隔为100m，地表逐点激发；
[0088]
为了模拟真实地震数据含有噪声和缺失道的影响，在理论地震数据集d_sei的基础上随机对15％的地震数据加入强度不等的高斯噪声，噪声标准差系数σ最大为0.5，另对15％的地震数据进行随机缺失地震道处理，最大随机缺失比为35％，这样得到观测地震数据d_sei_obs，扩充数据集的同时可以提高网络模型的鲁棒性；
[0089]
步骤二：先验知识处理和引入
[0090]
首先对于震源位置和检波器位置信息采用one-hot编码的方式引入到网络中；由于震源和检波器位置信息是地震数据和速度模型关联的特征，首先设置两行编码向量，一行表示震源位置信息，一行表示检波器位置信息，将其拼接到地震数据集d_sei_obs中，记为数据集d_sei_obs_pos；然后设置大小和速度模型一致的两个二维矩阵，通过one-hot编码方式一个反映震源位置，一个反映检波器位置,得到数据集d_vel_pos；其次测井数据作为地震勘探中可靠的先验知识，可以作为反演方法的地质先验信息约束。通过两口竖井的测井资料获取到速度模型数据集d_vel中第30列和第70列，并使用最近邻插值算法和高斯模糊操作生成粗糙的地下初始速度模型作为初始速度先验知识,记为数据集d_init；
[0091]
步骤三：网络输入与标签设定
[0092]
以步骤二中得到的加入震源和检波器位置信息的地震数据d_sei_obs_pos作为第一个编码器的输入，反映震源和检波器位置的数据集d_vel_pos和生成粗糙的地下初始速度模型数据集d_init作为第二个编码器的输入，以对应的速度模型数据集d_vel作为标签，为训练神经网络模型g提供数据；
[0093]
步骤四：基于物理先验约束的速度反演网络g构建
[0094]
网络模型g分为基于数据驱动的速度反演子网和基于rnn正演约束子网组成；基于数据驱动的速度反演子网由一个双编码器和一个解码器组成；基于rnn正演约束子网结合全波形反演理论和rnn网络，实现了物理规律与神经网络的有效结合；
[0095]
网络模型设定为
[0096]
(1)、数据驱动的速度反演子网gd[0097]
该子网的两个编码器分别输入加入震源和检波器位置信息的地震数据和生成粗糙的地下初始速度模型；通过测井先验提供大尺度波速模型，而地震数据的波形信息通过另一组编码器网络提供小尺度特征信息；通过两个编码网络将地震数据特征压缩至一维向量形式，格式为[bs，c，1，1]，其中bs是批处理大小、c为特征通道数，然后将提取的特征输入解码器得到网络的反演结果；
[0098]
(2)、rnn正演约束子网g
p
[0099]
将基于数据驱动的速度反演子网得到的反演结果作为rnn正演约束子网的输入，结合全波形反演理论和rnn网络处理时序数据的特征，将以波动方程驱动的波场延拓规律与rnn网络结合，通过卷积方式计算波场拉普拉斯值并利用隐藏状态层保存上个时刻的波场信息u,实现地震正演的过程，并结合神经网络自动微分算法和adam最优化算法实现自动求取梯度并更新参数模型；
[0100]
步骤五：损失函数设定
2080，深度学习网络模型使用pytorch0.4和python2.7搭建，波源子变换工具箱运行环境采用matlab2017b。实施的具体过程如下：
[0118]
1、地震数据训练集预处理
[0119]
1.1、读取mat文件
[0120]
通过工具包scipy.io的loadmat函数读取理论地震数据集d_sei和速度模型数据集d_vel中mat格式文件。调用方式为：
[0121]
dx＝loadmat(os.path.join(directory,filename_x,f"a.mat"))['f']
[0122]
dy＝loadmat(os.path.join(directory,filename_y,f"b.mat"))['f']
[0123]
将读取的地震数据存储在x，读取的速度模型存储在y，其中a.mat为地震数据文件名称，b.mat为速度模型文件名称，f为文件在matlab里的变量名称。
[0124]
1.2、地震数据的预处理
[0125]
为了模拟真实地震数据含有噪声和缺失道的影响，在理论地震数据集d_sei的基础上随机对15％的地震数据加入强度不等的高斯噪声，噪声标准差系数σ最大为0.5，另对15％的地震数据进行随机缺失地震道处理，最大随机缺失比为35％，得到观测地震数据d_sei_obs，关键代码如下：
[0126]
noise＝np.random.normal(0,(rand/30)*dx1.max(),(dx1.shape[0],1024,100))
[0127]
dx1＝dx1+noise
[0128]
rand1＝np.random.randint(100,size＝35)
[0129]
dx1[:,:,rand1]＝0.0
[0130]
然后使用one-hot编码将震源和检波器位置信息拼接到观测地震数据d_sei_obs中得到d_sei_obs_pos。
[0131]
1.3、测井数据制作
[0132]
设置大小和速度模型一致的两个二维矩阵，通过one-hot编码方式一个反映震源位置，一个反映检波器位置,得到数据集d_vel_pos；其次，通过两口竖井的测井资料获取到真实地震速度模型中第30列和第70列，并使用最近邻插值算法和高斯模糊操作生成粗糙的地下初始速度模型作为初始速度先验知识,记为数据集d_init。
[0133]
2、样本标签准备：
[0134]
以预处理得到的加入震源和检波器位置信息的地震数据d_sei_obs_pos作为第一个编码器的输入，反映震源和检波器位置的数据集d_vel_pos和生成粗糙的地下初始速度模型数据集d_init作为第二个编码器的输入，以对应的速度模型数据集d_vel作为标签，为训练神经网络模型g提供数据；
[0135]
3、设计基于物理约束的神经网络模型g的网络结构：本发明的网络模型可以划分为两部分，分别为速度反演子网和rnn正演约束子网。
[0136]
3.1、速度反演子网gd设计
[0137]
该子网络结合注意力机制和u-net网络，加强对先验知识的利用，主要包括卷积、批归一化函数、下采样、上采样和激活函数等操作，涉及到pytorch中相关的标准函数说明如下：
[0138]
构造卷积层时使用pytorch中定义的标准函数：nn.conv2d(in_channels,out_
channels,kernel_size,stride,padding)，其中in_channels表示输入通道数，out_channels表示输出通道数，stride指定了卷积核滑动的步长，padding为边缘填充的尺寸。
[0139]
批归一化函数batchnorm2d(num_features)，不仅可以加快模型的收敛速度，也能使数据在进行leakyrelu之前不会因为数据过大而导致网络性能的不稳定，参数num_features为其中特征的数量；
[0140]
leakyrelu层函数定义：nn.leakyrelu()，将前卷积层的输出结果输入到leakyrelu激活函数，调用函数以执行非线性映射。
[0141]
为了简化实施过程，首先将卷积、批归一化和激活函数3个操作定义为一组，每两组作为一个卷积块，每个编码器是由6个卷积块组成，将地震数据特征压缩至一维向量形式。同时，解码器也是由6个卷积块组成，将提取的特征解码至相应大小。卷积块的关键代码如下：
[0142][0143]
同时，利用attention_block模块，将编码器和解码器的对应层连接，attention_block模块关键代码如下：
[0144][0145]
3.2、rnn正演约束子网g
p
设计
[0146]
将速度反演子网的输出，作为rnn正演约束子网的输入。正演模拟过程就可以看作是对波场拉普拉斯值和波场值间的简单计算，对波场拉普拉斯值的计算可通过神经网络中的卷积操作实现，卷积操作的卷积核通常由有限差分格式和速度模型维度确定。对于二维波速模型的空间二阶有限差分格式，其卷积核为矩阵每层网络以前两个时刻的波场为输入，并引入当前时刻的震源函数值，经过计算得到新的时刻的波场值以及对应的观测地震数据。循环神经网络的关键代码为：h1，h2＝self.cell(h1，h2，c，s1)；其中，h1和h2为前两个时刻的波场值，c为速度值，s1为当前时刻的震源值。
[0147]
3.3、整体网络结构定义
[0148]
将上述基于数据驱动的速度反演子网和rnn正演约束网络串联起来，实现基于物理规律和先验知识约束的地震速度反演网络。该网络将反演网络输出的预测速度模型作为rnn正演约束网络的可训练参数输入使两个网络串联起来，在训练中经过rnn正演约束网络
得到的模拟地震数据和观测地震数据计算损失，实现对速度反演子网的梯度计算和参数更新，最终训练结束后，以反演子网的输出预测结果作为整体网络的最终结果。
[0149]
4、设计损失函数：
[0150]
整体网络的损失函数为：
[0151]
l＝l(x，y)+l(d
true
，d
pre
)
[0152]
x为速度模型数据集d_vel，y为速度反演子网的输出，d
true
和d
pre
分别为理论地震数据d_sei和rnn正演约束子网得到的预测地震数据。
[0153]
4.1、速度反演子网损失函数
[0154]
l(x，y)由均方误差函数和结构相似性函数(structural similarity，ssim)共同衡量x和y的差异，具体形式为：
[0155][0156][0157]
其中，μ
x
、μy表示x和y内样本元素的平均值，表示其方差，σ
xy
表示其协方差，c1和c2表示计算稳定性常数，n表示批样本数量。关键代码如下：
[0158]
loss_ssim＝content_criterion[0](sures_batch，hires_batch)
[0159]
loss_mse＝content_criterion[1](sures_batch，hires_batch)
[0160]
loss_g＝loss_mse-10*loss_ssim
[0161]
4.2、rnn正演约束子网损失函数
[0162]
l(d
true
，d
pre
)采用mse损失函数衡量模拟地震数据和观测地震数据的误差，结果为：
[0163][0164]
关键代码如下：loss_r＝crite_r(lores_batch，seismic_data)
[0165]
5、训练并保存网络模型：
[0166]
考虑到rnn网络内存消耗较大，且训练两个串联网络调参难度较大，本文基于迁移学习的思想，采用预训练的先验反演子网，结合rnn正演约束网络进行联合训练。整体训练流程为，先使用反演子网在数据集上训练迭代70次，然后将预训练的模型结合rnn正演约束网络继续训练迭代130次。最终训练结束后，以反演子网的输出预测结果作为整体网络的最终结果。
[0167]
6、测试网络模型性能：
[0168]
网络训练过程中的ssim和mse如图3和图4所示。测试在不同缺失道和和不同等级的噪声干扰下，有物理约束的网络模型a和无物理约束的网络模型b，对断层模型和盐丘模型的反演效果。同时，为说明本方法的有效性，将基于物理规律约束的速度反演网络和u-net网络反演方法以及inversionnet反演网络进行对比。
[0169]
实施效果：
[0170]
图5和图6给出模型a和模型b在测试集上不同程度缺失道测试结果对比，图5为不
同缺失道断层对比，其中第一行为无物理约束的网络，第二行为含有物理约束的网络；图6为不同缺失道盐丘对比，其中第一行为无物理约束的网络，第二行为含有物理约束的网络。可以看出针对缺失采样测试数据，模型a在缺失比例超过20％时保持了相对较好的结果，而模型b则出现细节上的变化，当缺失比例达到最大时有物理约束的网络虽有一定的差异产生但仍优于没有物理约束的网络。
[0171]
图7和图8给出模型a和模型b在测试集上不同噪声等级测试结果对比，图7为不同噪声等级断层对比，其中第一行为无物理约束的网络，第二行为含有物理约束的网络；图8为不同噪声等级盐丘对比，其中第一行为无物理约束的网络，第二行为含有物理约束的网络。可以看出在不同噪声数据测试中同样表明模型a对地震数据的噪声干扰有较好的鲁棒性。
[0172]
三种方法在测试集上表现如图9所示，第一行为标签数据；第二行为本文算法结果；第三行为u-net算法结果；第四行为inversion算法结果，从图中可以看出三种网络在简单层状模型和盐丘模型上重建较好，但在断层模型上，inversionnet网络在深层和断层出较为模糊，断层结构不明显；而u-net在断层细节处恢复不理想。给出三种网络模型在断层速度模型上第70列的速度值对比结果如图10所示，可见发明和标签的拟合度最好，以上结果均说明本方法的有效性。

技术特征：
1.一种基于物理约束神经网络的地震速度反演方法，其特征在于：地震速度反演方法包括以下步骤：步骤一、训练数据集准备速度反演任务建立由地震数据到速度模型的非线性映射建模，需要大量成对的地震数据-速度模型样本，但真实速度模型往往很难获得；采用人工生成二维速度模型(简单层状模型、断层模型和盐丘模型等)的方法扩充数据得到速度模型数据集d_vel；人工生成速度模型过程中，通过正弦、余弦函数以及线性函数的叠加生成层数在layer_m～layer_n之间，且层速度在velocity_m～velocity_n范围内的简单层状速度模型，速度模型尺寸为velocity_len
×
velocity_wid，横向和纵向的网格大小均为grid_size；在此基础上生成断层长度、位置随机及倾斜角度在angle_m～angle_n的断层模型；同时将大小为salt_dome_size个网格点、形状随机、速度为salt_dome_vel的盐丘随机加入断层模型中得到盐丘体模型；获得足够的速度模型后，采用有限差分正演的方法生成对应的理论地震数据集d_sei；对于震源和观测系统设置，采用地表放炮和地表观测方式，检波器在地表水平间隔每detector_interval_size放置一个(即地表每个网格点均有)，采样时间间隔为t，总采集时间为t，即接收到的地震道集记录有t/t个数据；震源子波采用主频为frequency_size的雷克子波，共shot_size炮，震源水平间隔为source_interval_size，地表逐点激发；为了模拟真实地震数据含有噪声和缺失道的影响，在理论地震数据集d_sei的基础上随机对dataset_percent_a的地震数据加入强度不等的高斯噪声，噪声标准差系数σ最大为noise_dev，另对dataset_percent_b的地震数据进行随机缺失地震道处理，最大随机缺失比为dele_percent，这样得到观测地震数据d_sei_obs，扩充数据集的同时可以提高网络模型的鲁棒性；步骤二：先验知识处理和引入对于震源位置和检波器位置信息采用one-hot编码的方式引入到网络中；由于震源和检波器位置信息是地震数据和速度模型关联的特征，设置两行编码向量，一行表示震源位置信息，一行表示检波器位置信息，将其拼接到地震数据集d_sei_obs中，记为数据集d_sei_obs_pos；然后设置大小和速度模型一致的两个二维矩阵，通过one-hot编码方式一个反映震源位置，一个反映检波器位置,得到数据集d_vel_pos；测井数据作为地震勘探中可靠的先验知识，可以作为反演方法的地质先验信息约束；通过well_num口竖井的测井资料获取到速度模型数据集d_vel中第col1列到第col
well_num
列，并使用最近邻插值算法和高斯模糊操作生成粗糙的地下初始速度模型作为初始速度先验知识,记为数据集d_init；步骤三：网络输入与标签设定以步骤二中得到的加入震源和检波器位置信息的地震数据d_sei_obs_pos作为第一个编码器的输入，反映震源和检波器位置的数据集d_vel_pos和生成粗糙的地下初始速度模型数据集d_init作为第二个编码器的输入，以对应的速度模型数据集d_vel作为标签，为训练神经网络模型g提供数据；步骤四：基于物理先验约束的速度反演网络g构建网络模型g分为基于数据驱动的速度反演子网和基于rnn正演约束子网组成；基于数据驱动的速度反演子网由一个双编码器和一个解码器组成；基于rnn正演约束子网结合全波
形反演理论和rnn网络，实现了物理规律与神经网络的有效结合；网络模型设定为(1)、数据驱动的速度反演子网g
d
速度反演子网g
d
的两个编码器分别输入加入震源和检波器位置信息的地震数据和生成粗糙的地下初始速度模型；通过测井先验提供大尺度波速模型，而地震数据的波形信息通过另一组编码器网络提供小尺度特征信息；通过两个编码网络将地震数据特征压缩至一维向量形式，格式为[bs，c，1，1]，其中bs是批处理大小、c为特征通道数，然后将提取的特征输入解码器得到网络的反演结果；(2)、rnn正演约束子网g
p
将基于数据驱动的速度反演子网得到的反演结果作为rnn正演约束子网的输入，结合全波形反演理论和rnn网络处理时序数据的特征，将以波动方程驱动的波场延拓规律与rnn网络结合，通过卷积方式计算波场拉普拉斯值并利用隐藏状态层保存上个时刻的波场信息u,实现地震正演的过程，并结合神经网络自动微分算法和adam最优化算法实现自动求取梯度并更新参数模型；步骤五：损失函数设定整体网络的损失函数为：l＝l(x，y)+l(d
true
，d
pre
)x为速度模型数据集d_vel，y为速度反演子网的输出，d
true
和d
pre
分别为理论地震数据d_sei和rnn正演约束子网得到的预测地震数据；其中l(x，y)由均方误差函数和结构相似性函数(structural similarity，ssim)共同衡量x和y的差异,具体形式为：衡量x和y的差异,具体形式为：其中,μ
x
、μ
y
表示x和y内样本元素的平均值，表示其方差，σ
xy
表示其协方差，c1和c2表示计算稳定性常数，n表示批样本数量；ssim是测量两个样本之间相似性的指标，主要关注边缘和纹理相似性，ssim通常取值范围为[0,1]，为1时表示两个样本完全相似，因此其取值范围从0开始越接近1越好，为了满足神经网络最小化损失函数的优化目标，采用负数形式满足要求；l(d
true
，d
pre
)采用mse损失函数衡量模拟地震数据和观测地震数据的误差,结果为：步骤六：预训练数据驱动的速度反演子网考虑到rnn网络内存消耗较大，且训练两个串联网络调参难度较大，基于迁移学习的思想，使用上述步骤二得到的加入先验知识的数据集预训练速度反演子网，将预训练好的网络模型保存；
步骤七：结合rnn正演约束网络联合训练将预训练的速度反演子网和rnn正演约束子网联合起来训练迭代，实现基于物理规律和先验知识约束的地震速度反演网络；该网络将反演网络输出的预测速度模型作为rnn正演约束网络的可训练参数输入使两个网络联通起来，在训练中经过rnn正演约束网络得到的预测地震数据和理论地震数据计算损失，实现对速度反演子网的梯度计算和参数更新，训练结束后，以反演子网的输出预测结果作为整体网络的最终速度反演结果。

技术总结
本发明属于地球科学技术与人工智能交叉领域，具体涉及一种基于物理约束神经网络的地震速度反演方法，该方法的步骤为1、训练数据集准备；2、先验知识处理和引入；3、网络输入与标签设定；4、基于物理先验约束的速度反演网络构建；5、损失函数设定；6、预训练数据驱动的速度反演子网；7、结合RNN正演约束网络联合训练。将数据驱动的神经网络和物理规律驱动的RNN正演约束网络结合起来实现基于数据和模型双驱动的速度反演网络，并结合迁移学习的思想先预训练数据驱动的速度反演网络，然后结合RNN正演约束网络进行联合训练，最终网络在较少样本的情况下实现高精度的速度反演。情况下实现高精度的速度反演。情况下实现高精度的速度反演。


技术研发人员：张岩 孟德聪 周一帆 张一鸣 张弛 宋利伟
受保护的技术使用者：东北石油大学
技术研发日：2023.06.26
技术公布日：2023/9/20