CN115032694A - 一种vsp初至旅行时层析成像方法及系统 - Google Patents

一种vsp初至旅行时层析成像方法及系统 Download PDF

Info

Publication number
CN115032694A
CN115032694A CN202210392912.2A CN202210392912A CN115032694A CN 115032694 A CN115032694 A CN 115032694A CN 202210392912 A CN202210392912 A CN 202210392912A CN 115032694 A CN115032694 A CN 115032694A
Authority
CN
China
Prior art keywords
vsp
seismic
data
equation
model
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.)
Pending
Application number
CN202210392912.2A
Other languages
English (en)
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 University of Petroleum East China
Hainan Institute of Zhejiang University
Original Assignee
China University of Petroleum East China
Hainan Institute of Zhejiang 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 China University of Petroleum East China, Hainan Institute of Zhejiang University filed Critical China University of Petroleum East China
Priority to CN202210392912.2A priority Critical patent/CN115032694A/zh
Publication of CN115032694A publication Critical patent/CN115032694A/zh
Pending legal-status Critical Current

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/38—Seismology; Seismic or acoustic prospecting or detecting specially adapted for water-covered areas
    • 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/301—Analysis for determining seismic cross-sections or geostructures

Landscapes

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

Abstract

本发明涉及一种VSP初至旅行时层析成像方法及系统，属于海洋矿产资源勘探和海底地层探测技术领域。先对地面地震剖面进行处理，得到地面地震构造信息数据，以构建初始VSP速度模型。再根据初始VSP速度模型计算初至旅行时计算值，将初至旅行时测量值和初至旅行时计算值进行相减，得到时间残差矩阵，以时间残差矩阵作为输入，对反演方程进行求解，得到更新后VSP速度模型，直至更新后VSP速度模型满足预定精度，从而利用地面地震构造信息数据作为先验知识对VSP速度模型的反演过程进行约束，提高VSP初至旅行时层析成像的精度。

Description

一种VSP初至旅行时层析成像方法及系统
技术领域
本发明涉及海洋矿产资源勘探和海底地层探测技术领域，特别是涉及一种地面地震构造信息数据约束下的VSP初至旅行时层析成像方法及系统。
背景技术
地震勘探是指激发人工震源所引起的地震波在地下介质中传播，由于地下介质弹性和密度的差异出现透射、反射、折射等现象，将地下介质的层位、岩性信息携带到地面，地球物理学家对地震波信号进行处理和解释，分析地震波在地下介质中的传播规律并进一步推断地下岩层的性质和形态。地震勘探的方法众多，根据震源和检波点的分布位置可以分为：地面地震勘探方法、垂直地震剖面(VSP)方法、井间地震方法。地面地震是指陆上地面、海洋拖缆和海底节点地震，地面地震勘探方法的震源和检波点近似在同一水平方向(以下称为水平面)，震源布置在地面或者海面以下(拖缆)，检波点布置在地面、海面以下(拖缆)或者海底，而VSP方法是由水平面震源激发地震波，固定在井壁的检波器接收地震记录的勘探方法。与地面地震方法相比，VSP地震勘探方法有着更高的分辨率和信噪比，可以提供钻井附近地层的空间和时间信息，更加细致地描述井旁地层地震构造。
VSP地震数据的观测方式减少了旅行时和射线传播路线，降低了近地表低速带引起的能量衰减，地震记录质量较高，初至旅行时易拾取，层析成像技术是建立VSP速度模型的有效途径。然而，VSP地震数据的数据量相对于地面地震数据较少，检波点的范围较小，目标成像范围也仅限于井周围，故综合利用先验信息约束提高反演精度仍然是VSP地震数据反演研究的重点方向，也是生产实践中急需解决的难点问题。
基于此，亟需一种能够提高反演精度的VSP初至旅行时层析成像方法及系统。
发明内容
本发明的目的是提供一种VSP初至旅行时层析成像方法及系统，利用地面地震构造信息数据作为先验知识对VSP速度模型的反演过程进行约束，提高VSP初至旅行时层析成像的精度。
为实现上述目的，本发明提供了如下方案：
一种VSP初至旅行时层析成像方法，所述层析成像方法包括：
对地面地震剖面进行处理，得到地面地震构造信息数据，并根据所述地面地震构造信息数据构建各向异性扩散方程；
以一维速度数据作为输入，对所述各向异性扩散方程进行求解，得到初始VSP速度模型；所述一维速度数据为VSP速度数据或者测井速度数据；
提取VSP地震数据中的初至旅行时测量值；
根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值；
将所述初至旅行时测量值和所述初至旅行时计算值进行相减，得到时间残差矩阵；
以所述时间残差矩阵作为输入，对反演方程进行求解，得到更新后VSP速度模型；
判断所述更新后VSP速度模型是否满足预定精度；
若是，则以所述更新后VSP速度模型作为最终的VSP速度模型；
若否，则以所述更新后VSP速度模型作为下一迭代中的初始VSP速度模型，返回“根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值”的步骤。
一种VSP初至旅行时层析成像系统，所述层析成像系统包括：
各向异性扩散方程构建模块，用于对地面地震剖面进行处理，得到地面地震构造信息数据，并根据所述地面地震构造信息数据构建各向异性扩散方程；
初始VSP速度模型求解模块，用于以一维速度数据作为输入，对所述各向异性扩散方程进行求解，得到初始VSP速度模型；所述一维速度数据为VSP速度数据或者测井速度数据；
时间残差矩阵计算模块，用于提取VSP地震数据中的初至旅行时测量值；根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值；将所述初至旅行时测量值和所述初至旅行时计算值进行相减，得到时间残差矩阵；
迭代更新模块，用于以所述时间残差矩阵作为输入，对反演方程进行求解，得到更新后VSP速度模型；判断所述更新后VSP速度模型是否满足预定精度；若是，则以所述更新后VSP速度模型作为最终的VSP速度模型；若否，则以所述更新后VSP速度模型作为下一迭代中的初始VSP速度模型，返回“根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值”的步骤。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本实施例用于提供一种VSP初至旅行时层析成像方法及系统，先对地面地震剖面进行处理，得到地面地震构造信息数据，并根据地面地震构造信息数据构建各向异性扩散方程，然后以一维速度数据作为输入，对各向异性扩散方程进行求解，得到初始VSP速度模型。再提取VSP地震数据中的初至旅行时测量值，根据初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值，将初至旅行时测量值和初至旅行时计算值进行相减，得到时间残差矩阵，以时间残差矩阵作为输入，对反演方程进行求解，得到更新后VSP速度模型，直至更新后VSP速度模型满足预定精度，则以更新后VSP速度模型作为最终的VSP速度模型，从而利用地面地震构造信息数据作为先验知识对VSP速度模型的反演过程进行约束，提高VSP初至旅行时层析成像的精度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例1所提供的层析成像方法的方法流程图；
图2为本发明实施例1所提供的层析成像方法的方法原理图；
图3为本发明实施例1所提供的二维层状模型示意图；
图4为本发明实施例1所提供的二维层状模型的叠加记录示意图；
图5为本发明实施例1所提供的VSP地震记录及初至旅行时拾取示意图；
图6为本发明实施例1所提供的井筒位置提取的测井速度示意图，其中，图6(a)为来自模型的真实速度，图6(b)为加入随机噪声后所得的声波测井速度，图6(c)为从图6(b)中得到的平滑速度；
图7为本发明实施例1所提供的平滑测井速度沿构造方向扩散构建的初始模型示意图；
图8为本发明实施例1所提供的层析反演得到的速度模型示意图；
图9为本发明实施例1所提供的层析反演的速度更新量示意图；
图10为本发明实施例2所提供的层析成像系统的系统框图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种VSP初至旅行时层析成像方法及系统，利用地面地震构造信息数据作为先验知识对VSP速度模型的反演过程进行约束，提高VSP初至旅行时层析成像的精度。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
实施例1：
经过多年发展，VSP技术的应用研究不断深入，从零偏移距VSP技术发展到变偏移距VSP技术、三维VSP技术、多分量VSP技术、逆VSP技术以及Walkaway VSP技术等，且随着DAS技术的发展，VSP技术得到了更广泛的推广和应用。然而，VSP地震数据的数据量相对于地面地震较少，检波点的范围较小，目标成像范围也仅限于井周围，故综合利用先验信息约束提高反演精度仍然是VSP地震数据反演研究的重点方向，也是生产实践中急需解决的难点问题。
VSP地震数据具有高精度、高分辨率的优势，但是VSP勘探范围有限，很难准确反演介质速度模型，而地面地震数据的精度虽然不如VSP地震数据，但勘探范围广，故利用地面地震数据与VSP地震数据进行联合层析成像，可以实现各种数据间的优势互补，从而得到高质量的反演结果。目前关于地面地震数据和VSP地震数据的联合反演方法的研究包括：(1)利用测井、VSP和地面地震数据旅行时层析信息来反演地震速度模型；(2)将地面地震数据、VSP地震数据和井间地震记录偏移到同一成像空间，以进行立体地震成像，得到井眼附近区域的高质量成像；(3)利用VSP地震数据和地面地震数据联合层析反演得到地下介质的速度；(4)利用VSP地震数据的初至旅行时反演各向异性参数，辅助建立地面地震各向异性介质模型；(5)开展井间地震、VSP和地面地震数据的联合成像，对不同尺度地震数据进行叠前深度偏移，应用归一化处理后的自适应叠加方法实现不同尺度数据间的联合成像。众多联合反演方法试图解决VSP地震数据和地面地震数据各自的不足，通过联合地面地震数据来提高VSP地震数据建模的精度。本实施例仍然基于综合利用地面地震和VSP地震数据的思路，提出一种新的地面地震构造信息数据约束下的VSP初至旅行时层析成像方法，主要就利用各向异性扩散方程建模和层析成像预条件正则化进行深入研究，从而达到利用先验的地面地震构造信息数据提高VSP层析成像精度的目的。
本实施例用于提供一种地面地震构造信息数据约束下的VSP初至旅行时层析成像方法，通过引入先验的地面地震构造信息数据来提高VSP地震数据层析成像的精度，如图1和图2所示，该层析成像方法包括：
S1：对地面地震剖面进行处理，得到地面地震构造信息数据，并根据所述地面地震构造信息数据构建各向异性扩散方程；
本实施例的地面地震剖面包括时间域地面地震剖面和深度域地面地震剖面，则对地面地震剖面进行处理，得到地面地震构造信息数据可以包括：
(1)判断地面地震剖面是否为时间域地面地震剖面；
(2)若是，则先将时间域地面地震剖面转化为深度域地面地震剖面，再对深度域地面地震剖面进行倾角扫描，得到地面地震构造信息数据；
(3)若否，即地面地震剖面为深度域地面地震剖面，则直接对深度域地面地震剖面进行倾角扫描，得到地面地震构造信息数据。
本实施例对深度域地面地震剖面进行倾角扫描，得到地面地震构造信息数据，若地面地震剖面为时间域地面地震剖面则需要先转化为深度域地面地震剖面，地面地震构造信息数据是指地面地震剖面的构造信息，用来约束层析成像过程。地面地震构造信息数据使用标准的SEGY格式，标准SEGY格式的文件一般包括三部分，第一部分是EBCDIC文件头，用来保存一些对地震数据体进行描述的信息；第二部分是二进制文件头，用来存储描述SEGY文件的关键信息，如数据格式、采样点数、采样间隔等信息；第三部分是实际的地震道，每条地震道都包含道头信息和地震道数据。或者，地面地震构造信息数据使用包含纵测线(Xline)和横测线(Inline)位置坐标、纵向采样间隔和采样点信息的其他数据格式。
本实施例所用的VSP地震数据和地面地震数据这两种地震数据之间的关系是地面地震成像剖面测线包含或者邻近VSP井位，可以在VSP地震数据旅行时层析成像过程中利用地面地震剖面的构造信息来提高反演精度，因此平滑步骤不能破坏构造信息，即必须使用保护地面地震构造信息的平滑方法。基于偏微分方程的各向异性扩散方程是一类自适应的平滑方法，可以根据不同图像内容采取不同的平滑方式，方程的解就是平滑的结果。
最早的各向异性扩散方程可以追溯到Perona和Malik提出的PM模型，其形式如下：
Figure BDA0003596250910000061
其中，u(x，z；t)为(x，z)位置点处，t时刻时的地震图像；div为散度算子；g(·)为非负递减的扩散系数；
Figure BDA0003596250910000062
为梯度算子；u0(x,z)为原始地震图像。
Weicket通过对图像数据的结构进行分析，将扩散系数转换为与方向性特征有关的扩散张量，构建了如下张量型扩散方程的模型：
Figure BDA0003596250910000063
上式即为本实施例所用的各向异性扩散方程。其中，u(x，z；t)为由一维速度数据所得到的t时刻，(x，z)位置点处的地震图像；div为散度算子；D为扩散张量，也可称之为张量型扩散系数，它的元素是根据结构张量Sρ所提取的图像局部结构信息来设计的，用以控制扩散的方向以及在对应方向上的扩散速率；Sρ为结构张量；
Figure BDA0003596250910000064
为梯度算子；u0(x，z)为由一维速度数据所得到的(x，z)位置点处的初始地震图像。
本实施例即利用上述得到的地面地震构造信息数据来求解扩散张量，以构建完整的各向异性扩散方程，则根据地面地震构造信息数据构建各向异性扩散方程可以包括：
(1)计算地面地震构造信息数据的结构张量；
Figure BDA0003596250910000065
uσ＝u*Gσ，其中，Sρ为地面地震构造信息数据的结构张量；Gρ表示尺度为ρ的Gaussion函数；
Figure BDA0003596250910000066
为梯度算子；u为由地面地震构造信息数据所得到的(x，y)位置点处，t时刻的地震图像；Gσ表示尺度为σ的Gaussion函数，用来避免计算梯度时噪声的影响；Gρ与梯度平方张量
Figure BDA0003596250910000067
卷积的目的在于考虑周围信息，避免在确定边缘方向时由于方向相同、符号相反而相互抵消。
结构张量使得二维图像中的每个点都对应一个与一阶偏微分信息有关的半正定矩阵，则结构张量为：
Figure BDA0003596250910000068
其中，Sρ为结构张量；Gρ为尺度为ρ的Gaussion函数；
Figure BDA0003596250910000071
为梯度算子；uσ＝u*Gσ，u为由地面地震构造信息数据所得到的(x，y)位置点处的地震图像；Gσ为尺度为σ的Gaussion函数。
(2)根据矩阵特征值的分解定理，对结构张量进行特征值分解，得到结构张量的第一特征值和第一特征向量；
特征值分解后，可以将结构张量表示为：
Figure BDA0003596250910000072
其中，v1,v2为第一特征向量；λ1,λ2为第一特征向量对应的第一特征值；v1表示平行
Figure BDA0003596250910000073
的方向，即变化率最大的方向，变化率为
Figure BDA0003596250910000074
v2表示垂直于
Figure BDA0003596250910000075
的方向，即变化率最小的方向，变化率为
Figure BDA0003596250910000076
第一特征向量的具体解形式如下：
Figure BDA0003596250910000077
对应的特征值为：
Figure BDA0003596250910000078
(3)对第一特征值进行归一化和正则化处理，得到第二特征值；以第二特征值和第一特征向量作为扩散张量的特征值和特征向量，构建各向异性扩散方程中的扩散张量，得到各向异性扩散方程。
扩散张量的特征向量与结构张量相同，但是特征值不同，特征值需要根据图像处理的具体需求来计算。Weickert根据边缘增强和相干增强两种图像处理要求给出了扩散张量的特征值，特征值u1和u2分别独立控制在v1和v2方向上的扩散行为，此时选取的特征值为：
Figure BDA0003596250910000079
通过上述公式得到的扩散张量，利用结构张量的特征向量对图像的变化进行准确定向，同时可以根据图像处理的需求定义沿特征方向变化的大小(即特征值的大小)，不再只依赖于图像的局部结构张量或局部结构。
区别于图像处理，地震勘探问题需要满足以下两个需求：扩散张量看作是结构张量的函数满足无穷光滑，扩散张量的特征值存在正下界。根据两个需求并结合地震勘探的实际问题，本实施例通过以下两个步骤来求解扩散张量的特征值：
(3.1)对第一特征值进行L2范数归一化，得到归一化后特征值；
Figure BDA0003596250910000081
其中，A、B为归一化后特征值，λ＝(λ1，λ2)是特征值向量。
(3.2)根据扩散张量特征值的下界对归一化后特征值进行正则化，得到第二特征值。
考虑到扩散张量特征值存在正下界，则归一化后特征值通过下式正则化到[α,1]范围内：
Figure BDA0003596250910000082
根据扩散张量和结构张量之间的关系，利用上述新的第二特征值和结构张量的特征向量可以构建正则扩散张量。在得到正则扩散张量后，即可建立得到各向异性扩散方程。
VSP地震勘探的观测范围较小导致VSP数据具有不完备性，更造成了层析成像的不适定性，因此在层析成像中加入先验信息是非常必要的。通过构造各向异性扩散方程的扩散张量，将地面地震构造信息引入到VSP地震数据层析成像中，使用先验信息约束层析成像。
S2：以一维速度数据作为输入，对所述各向异性扩散方程进行求解，得到初始VSP速度模型；所述一维速度数据为VSP速度数据或者测井速度数据；
将VSP速度数据和测井速度数据统一称为一维速度数据，则一维速度数据采用文本格式，包含井名，位置和深度以及速度信息。对一维速度数据(VSP速度数据或者测井速度数据)进行检查，去除不合理的异常值(负值等)，得到合理的一维速度数据，该合理的一维速度数据即为本实施例S2中的一维速度数据。以一维速度数据作为输入，对各向异性扩散方程进行求解的具体求解方式为：构建一个空矩阵，根据一维速度数据的位置信息将一维速度数据填充到矩阵中对应的位置，填充后的矩阵作为各向异性扩散方程中的初始图像u0，然后求解方程进行扩散，得到初始VSP速度模型。
本实施例通过深度域地面地震剖面得到地面地震构造信息数据，利用该地面地震构造信息数据来构建各向异性扩散方程中的扩散张量，以构建得到各向异性扩散方程，然后基于一维速度数据求解各向异性扩散方程，将井筒处的速度扩散为二维或三维的初始VSP速度模型，在得到初始VSP速度模型后，将其作为输入进行VSP初至旅行时层析成像。基于一维速度数据，使用各向异性扩散方程沿地面地震大尺度构造方向扩散建立初始VSP速度模型，能够利用地面地震构造信息数据为层析成像提供良好的初始模型。
S3：提取VSP地震数据中的初至旅行时测量值；
本实施例的VSP地震数据使用标准的SEGY格式，或者该VSP地震数据使用包含每道数据的初至波旅行时、炮号、震源点和检波点的坐标、时间采样间隔和采样点数信息的数据格式，震源点、检波点的位置可以用来确定地震波激发和接收的位置，VSP地震数据的初至波旅行时，与快速扫描法求取的旅行时作差，得到旅行时误差，进一步对反演方程进行求解。
S4：根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值；
具体的，S4可以包括：根据初始VSP速度模型构建程函方程；利用快速扫描法对所述程函方程进行求解，得到初至旅行时计算值。
更为具体的，在“高频近似”的假设下(即地震波的波长远小于介质的不均匀性，地震波能量主要沿射线传播)，弹性波方程可以被简化为程函方程，程函方程给出了任意一点的旅行时梯度的大小与慢度(速度的倒数)之间的关系，二维空间中程函方程表达式为：
Figure BDA0003596250910000091
其中，
Figure BDA0003596250910000092
为梯度算子，T为旅行时，T(x)为x位置点的旅行时，c(x)为地震波在介质中x位置点处的传播速度。直达波的旅行时梯度与慢度满足程函方程的表达波，且在震源点处满足T(xs)＝0，xs为震源的位置，
快速扫描法基于因果关系将旅行时场传播的方向分为有限个组，对于每一组分别利用Gauss-Seidel迭代方程求解离散化后的逆风差分格式的离散方程，逆风差分格式的离散方程为：
Figure BDA0003596250910000093
其中，ti,j为网格点(i，j)处的旅行时；
Figure BDA0003596250910000094
i，j分别为x，z方向的网格点数，si,j为网格点(i，j)处的慢度(速度的倒数)，h为空间步长，且
Figure BDA0003596250910000101
其中，f代指()+中的数据或者表达式，若f大于0，则(f)+＝f，若f小于或等于0，则(f)+＝0。
在二维空间中，可以将扫描方向分成四组(右上，左上，左下和右下)，初至值信息沿着四个方向进行传播，通过四组不同顺序的Gauss-Seidel迭代并结合逆风差分格式求解地震波的传播旅行时。
快速扫描法计算地震初至波旅行时的步骤为：
(1)初始化：在震源点上赋初值为0，震源点周围的四个点上赋值为精确值以避免一阶数值误差，震源点周围的四个点是指震源点上下左右四个点，可以通过求解程函方程获得旅行时的精确值，除震源点和震源点之外的其它网格点的初值设为一个较大值，使其远大于所有计算点最终能算出的旅行时，这些点上的旅行时在后续计算时会被更新；
(2)迭代：在四个方向上，使用Gauss-Seidel型迭代法求解离散的非线性方程组，在每一个旅行时未被固定的网格点，计算当前旅行时并与原值比较，取较小值作为当前值；
(3)判断收敛：如果满足收敛条件||tnew-told||≤δ，计算结束；如果不满足，那么继续迭代直至收敛为止。利用逆风差分格式的离散方程来求解除震源点和上述周围四个点之外的每一个网格点的当前走时。
S5：将所述初至旅行时测量值和所述初至旅行时计算值进行相减，得到时间残差矩阵；以所述时间残差矩阵作为输入，对反演方程进行求解，得到更新后VSP速度模型；
层析成像的基本原理是先给定一个初始慢度模型S(速度模型的倒数)，根据射线路径追踪，构造射线路径矩阵L，通过快速扫描法求解程函方程得到地震波初至旅行时计算值，将初至旅行时测量值与初至旅行时计算值相减得到旅行时残差矩阵ΔT，从而形成反演方程：
LΔS＝ΔT；
其中，L为射线路径矩阵；ΔS为慢度模型，慢度为速度的倒数；ΔT为时间残差矩阵。
求解上述大型稀疏方程组就可以得到介质慢度，即得到更新后VSP速度模型。该方程是一个病态的大型系数线性代数方程组，可用代数重建技术(ART)、瞬时迭代重建技术(SIRT)、阻尼最小二乘法(LSQR)等方法进行求解。
S6：判断所述更新后VSP速度模型是否满足预定精度；若是，则以所述更新后VSP速度模型作为最终的VSP速度模型；若否，则以所述更新后VSP速度模型作为下一迭代中的初始VSP速度模型，返回“根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值”的步骤，通过不断迭代得到最终的反演速度。
作为一种可选的实施方式，在得到更新后VSP速度模型，并在判断更新后VSP速度模型是否满足预定精度之前，本实施例的层析成像方法还包括：判断迭代次数是否达到预定次数；若是，则利用各向异性扩散方程对反演方程中的慢度模型进行预条件平滑，得到新的反演方程，并以新的反演方程作为下一迭代中的反演方程，同时将迭代次数置为0；若否，则不改变反演方程，继续迭代。
具体的，在模型更新过程中，可以利用地面地震构造信息数据构建各向异性扩散方程，对慢度模型进行预条件平滑：
ΔS′＝QΔS；
上式是各向异性扩散方程的示意表示形式，各向异性扩散方程就是实现上述过程的滤波。其中，Q为预条件算子，QΔS即表示利用各向异性扩散方程对模型ΔS进行平滑扩散。为了提高反演效率，不必每一次更新都通过各向异性扩散方程进行平滑，可以间歇地沿地面地震构造方向平滑模型，然后进入下一次模型更新。
本实施例利用预条件正则化沿地震构造方向对模型进行平滑扩散，能够引入地面地震构造信息来提高层析成像的质量，同时在迭代过程中间歇地对速度模型沿地震构造方向进行扩散平滑，以更高效率地构建高精度VSP速度模型。
本实施例先对深度域地面地震剖面进行倾角扫描得到地面地震构造信息，基于一维VSP速度数据或测井速度数据，使用各向异性扩散方程沿地面地震大尺度构造方向构建初始速度模型，使用快速扫描法求解程函方程得到地震初至旅行时，求解大型稀疏矩阵线性代数方程组反演速度模型，实现高效准确的VSP资料层析成像，利用地面地震构造信息进行预条件层析成像，在层析反演的迭代步骤中，间歇地对速度模型沿地面地震构造方向进行扩散平滑，以提高VSP资料层析成像的精度，从而利用地震先验构造信息提高VSP初至旅行时层析成像的精度。
以下，使用二维层状模型测试预条件VSP资料初至波旅行时层析方法的有效性，层状模型如图3所示，大小为9000m×6000m，密度为常量，速度包含多层起伏变化的界面约束的层位，模型速度从上到下分别为2500米/秒、3000米/秒、3400米/秒、4100米/秒和4500米/秒，模型构造复杂，地层起伏变化较大，对其进行层析成像也具有一定难度。VSP地震记录数据是通过速度模型(图3)正演模拟获取，共121炮，炮间距为50米，震源埋深为0.5米，0.6千米到6.6千米范围内均匀分布，井筒位于水平方向0.5千米处，接收点在1千米到5千米的深度范围内均匀分布。
图4是模拟的叠加剖面，主要用来生成速度模型所在区域的地下构造信息，该地面地震剖面是深度域地震剖面，实际生产中地震剖面大多是叠前时间偏移剖面，可以通过时深转换得到对应的深度域剖面。图5是基于速度模型获取的VSP地震记录，为了进行VSP地震数据的初至波旅行时层析成像，需要拾取初至旅行时，图中圆圈表示拾取的初至波到达时间，作为初至旅行时层析的输入数据。
为了模拟真实处理的步骤，图6表示在VSP井位置，图6(a)是来自模型的真实速度，图6(b)是加入随机噪声后所得的声波测井速度，图6(c)是从图6(b)中得到的平滑速度。本实施例基于一维速度模型利用各向异性扩散方程得到层析成像的初始模型，如图7所示，初始速度模型由于使用了地震叠加剖面的构造作为约束扩散而成，明显包含地面构造信息，整体上与真实速度场较接近。图8是基于初始模型进行层析反演得到的速度模型，图9是层析反演的更新量。层析成像更新结果与初始模型相比变化较明显，主要是构造边界处更新量较大，更新量中有叠加记录的构造痕迹。
实施例2：
本实施例用于提供一种VSP初至旅行时层析成像系统，如图10所示，所述层析成像系统包括：
各向异性扩散方程构建模块M1，用于对地面地震剖面进行处理，得到地面地震构造信息数据，并根据所述地面地震构造信息数据构建各向异性扩散方程；
初始VSP速度模型求解模块M2，用于以一维速度数据作为输入，对所述各向异性扩散方程进行求解，得到初始VSP速度模型；所述一维速度数据为VSP速度数据或者测井速度数据；
时间残差矩阵计算模块M3，用于提取VSP地震数据中的初至旅行时测量值；根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值；将所述初至旅行时测量值和所述初至旅行时计算值进行相减，得到时间残差矩阵；
迭代更新模块M4，用于以所述时间残差矩阵作为输入，对反演方程进行求解，得到更新后VSP速度模型；判断所述更新后VSP速度模型是否满足预定精度；若是，则以所述更新后VSP速度模型作为最终的VSP速度模型；若否，则以所述更新后VSP速度模型作为下一迭代中的初始VSP速度模型，返回“根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值”的步骤。
本说明书中每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种VSP初至旅行时层析成像方法，其特征在于，所述层析成像方法包括：
对地面地震剖面进行处理，得到地面地震构造信息数据，并根据所述地面地震构造信息数据构建各向异性扩散方程；
以一维速度数据作为输入，对所述各向异性扩散方程进行求解，得到初始VSP速度模型；所述一维速度数据为VSP速度数据或者测井速度数据；
提取VSP地震数据中的初至旅行时测量值；
根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值；
将所述初至旅行时测量值和所述初至旅行时计算值进行相减，得到时间残差矩阵；
以所述时间残差矩阵作为输入，对反演方程进行求解，得到更新后VSP速度模型；
判断所述更新后VSP速度模型是否满足预定精度；
若是，则以所述更新后VSP速度模型作为最终的VSP速度模型；
若否，则以所述更新后VSP速度模型作为下一迭代中的初始VSP速度模型，返回“根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值”的步骤。
2.根据权利要求1所述的层析成像方法，其特征在于，所述对地面地震剖面进行处理，得到地面地震构造信息数据具体包括：
判断所述地面地震剖面是否为时间域地面地震剖面；
若是，则将所述时间域地面地震剖面转化为深度域地面地震剖面；
对所述深度域地面地震剖面进行倾角扫描，得到地面地震构造信息数据。
3.根据权利要求1所述的层析成像方法，其特征在于，所述根据所述地面地震构造信息数据构建各向异性扩散方程具体包括：
计算所述地面地震构造信息数据的结构张量；
对所述结构张量进行特征值分解，得到所述结构张量的第一特征值和第一特征向量；
对所述第一特征值进行归一化和正则化处理，得到第二特征值；
以所述第二特征值和所述第一特征向量作为扩散张量的特征值和特征向量，构建各向异性扩散方程中的扩散张量，得到各向异性扩散方程。
4.根据权利要求3所述的层析成像方法，其特征在于，所述结构张量为：
Figure FDA0003596250900000021
其中，Sρ为结构张量；Gρ为尺度为ρ的Gaussion函数；
Figure FDA0003596250900000022
为梯度算子；uσ＝u*Gσ，u为由地面地震构造信息数据所得到的(x，y)位置点处的地震图像；Gσ为尺度为σ的Gaussion函数。
5.根据权利要求3所述的层析成像方法，其特征在于，所述对所述第一特征值进行归一化和正则化处理，得到第二特征值具体包括：
对所述第一特征值进行L2范数归一化，得到归一化后特征值；
根据扩散张量特征值的下界对所述归一化后特征值进行正则化，得到第二特征值。
6.根据权利要求1所述的层析成像方法，其特征在于，所述各向异性扩散方程为：
Figure FDA0003596250900000023
其中，u(x，z；t)为由一维速度数据所得到的t时刻，(x，z)位置点处的地震图像；div为散度算子，D为扩散张量；Sρ为结构张量；
Figure FDA0003596250900000024
为梯度算子；u0(x，z)为由一维速度数据所得到的(x，z)位置点处的初始地震图像。
7.根据权利要求1所述的层析成像方法，其特征在于，所述根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值具体包括：
根据所述初始VSP速度模型构建程函方程；
利用快速扫描法对所述程函方程进行求解，得到初至旅行时计算值。
8.根据权利要求1所述的层析成像方法，其特征在于，所述反演方程为：
LΔS＝ΔT；
其中，L为射线路径矩阵；ΔS为慢度模型，慢度为速度的倒数；ΔT为时间残差矩阵。
9.根据权利要求1所述的层析成像方法，其特征在于，在得到更新后VSP速度模型，并在判断所述更新后VSP速度模型是否满足预定精度之前，所述层析成像方法还包括：
判断迭代次数是否达到预定次数；
若是，则利用所述各向异性扩散方程对所述反演方程中的慢度模型进行预条件平滑，得到新的反演方程，并以所述新的反演方程作为下一迭代中的反演方程；同时将迭代次数置为0。
10.一种VSP初至旅行时层析成像系统，其特征在于，所述层析成像系统包括：
各向异性扩散方程构建模块，用于对地面地震剖面进行处理，得到地面地震构造信息数据，并根据所述地面地震构造信息数据构建各向异性扩散方程；
初始VSP速度模型求解模块，用于以一维速度数据作为输入，对所述各向异性扩散方程进行求解，得到初始VSP速度模型；所述一维速度数据为VSP速度数据或者测井速度数据；
时间残差矩阵计算模块，用于提取VSP地震数据中的初至旅行时测量值；根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值；将所述初至旅行时测量值和所述初至旅行时计算值进行相减，得到时间残差矩阵；
迭代更新模块，用于以所述时间残差矩阵作为输入，对反演方程进行求解，得到更新后VSP速度模型；判断所述更新后VSP速度模型是否满足预定精度；若是，则以所述更新后VSP速度模型作为最终的VSP速度模型；若否，则以所述更新后VSP速度模型作为下一迭代中的初始VSP速度模型，返回“根据所述初始VSP速度模型，利用快速扫描法计算得到初至旅行时计算值”的步骤。
CN202210392912.2A 2022-04-14 2022-04-14 一种vsp初至旅行时层析成像方法及系统 Pending CN115032694A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210392912.2A CN115032694A (zh) 2022-04-14 2022-04-14 一种vsp初至旅行时层析成像方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210392912.2A CN115032694A (zh) 2022-04-14 2022-04-14 一种vsp初至旅行时层析成像方法及系统

Publications (1)

Publication Number Publication Date
CN115032694A true CN115032694A (zh) 2022-09-09

Family

ID=83119877

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210392912.2A Pending CN115032694A (zh) 2022-04-14 2022-04-14 一种vsp初至旅行时层析成像方法及系统

Country Status (1)

Country Link
CN (1) CN115032694A (zh)

Similar Documents

Publication Publication Date Title
RU2693495C1 (ru) Полная инверсия волнового поля с компенсацией показателя качества
Wu et al. A procedure for 3D simulation of seismic wave propagation considering source‐path‐site effects: Theory, verification and application
US8363509B2 (en) Method for building velocity models for pre-stack depth migration via the simultaneous joint inversion of seismic, gravity and magnetotelluric data
CN105277978B (zh) 一种确定近地表速度模型的方法及装置
CN108802813B (zh) 一种多分量地震资料偏移成像方法及系统
US10234582B2 (en) Joint inversion of seismic data
US20110090760A1 (en) Full-waveform inversion in the traveltime domain
CN109884710B (zh) 针对激发井深设计的微测井层析成像方法
US20170335675A1 (en) Method To Predict Pore Pressure And Seal Integrity Using Full Wavefield Inversion
US11609349B2 (en) Determining properties of a subterranean formation using an acoustic wave equation with a reflectivity parameterization
CN106526674A (zh) 一种三维全波形反演能量加权梯度预处理方法
CN103513277B (zh) 一种地震地层裂隙裂缝密度反演方法及系统
CN109557582B (zh) 一种二维多分量地震资料偏移成像方法及系统
CN113740901B (zh) 基于复杂起伏地表的陆上地震数据全波形反演方法及装置
US11733413B2 (en) Method and system for super resolution least-squares reverse time migration
Golubev et al. Simulation of seismic wave propagation in a multicomponent oil deposit model
US20160161619A1 (en) Systems and methods for aligning a monitor seismic survey with a baseline seismic survey
WO2016001697A1 (en) Systems and methods for geologic surface reconstruction using implicit functions
Huang et al. 2D/3D seismic simultaneous inversion for the velocity and interface geometry using multiple classes of arrivals
CN104199088B (zh) 一种提取入射角道集的方法及系统
US6324478B1 (en) Second-and higher-order traveltimes for seismic imaging
Guo et al. Topography-dependent eikonal tomography based on the fast-sweeping scheme and the adjoint-state technique
Waheed et al. A holistic approach to computing first-arrival traveltimes using neural networks
CN115032694A (zh) 一种vsp初至旅行时层析成像方法及系统
CN113267830A (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