CN117687083A - 全波形反演方法、设备及存储介质 - Google Patents
全波形反演方法、设备及存储介质 Download PDFInfo
- Publication number
- CN117687083A CN117687083A CN202211105068.7A CN202211105068A CN117687083A CN 117687083 A CN117687083 A CN 117687083A CN 202211105068 A CN202211105068 A CN 202211105068A CN 117687083 A CN117687083 A CN 117687083A
- Authority
- CN
- China
- Prior art keywords
- seismic data
- seismic
- data
- matrix
- representing
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 51
- 239000011159 matrix material Substances 0.000 claims abstract description 83
- 208000037516 chromosome inversion disease Diseases 0.000 claims description 81
- 230000006870 function Effects 0.000 claims description 44
- 239000013598 vector Substances 0.000 claims description 31
- 238000006073 displacement reaction Methods 0.000 claims description 12
- 230000000875 corresponding effect Effects 0.000 claims description 10
- 230000002596 correlated effect Effects 0.000 claims description 4
- 238000013507 mapping Methods 0.000 claims description 4
- 238000010606 normalization Methods 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 14
- 230000005540 biological transmission Effects 0.000 description 11
- 238000012545 processing Methods 0.000 description 5
- 238000004590 computer program Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 230000008878 coupling Effects 0.000 description 2
- 238000010168 coupling process Methods 0.000 description 2
- 238000005859 coupling reaction Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000005314 correlation function Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
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
-
- 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
-
- 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供一种全波形反演方法、设备及存储介质，该方法包括：根据观测地震数据和合成地震数据构建匹配代价矩阵，匹配代价矩阵中的每一个矩阵元素用于表示合成地震数据中一个地震道数据与观测地震数据中一个地震道数据之间的相似性；根据匹配代价矩阵确定用于反映观测地震数据和合成地震数据误差的目标函数，根据目标函数得到伴随场震源；根据伴随场震源与正传波场确定梯度；根据梯度和预设迭代步长对初始速度模型进行迭代更新，直至满足预设收敛准则。将合成地震数据的地震道与观测地震数据的地震道进行匹配，充分考虑了地震道的空间特性，降低了对于初始速度模型的要求，不仅能够避免陷入局部极小值，而且能够提高反演精度。
Description
技术领域
本发明实施例涉及油气勘探技术领域，具体涉及一种全波形反演方法、设备及存储介质。
背景技术
工业勘探技术日趋成熟，从地震数据采集到最后油气藏层位识别已经形成了一整套完整的流程：首先进行地震数据采集，在地震数据采集完毕之后，对其进行数据预处理，去除噪声并标记初至波；然后建立初始速度模型，采用全波形反演方法从初始速度模型出发，重建高精度高保真的地下介质速度模型；最后基于重建的地下介质速度模型进行地震成像和地质解释。在整个地震勘探流程中，全波形反演方法的效果起着关键作用，直接决定了最后是否能成功识别和判断出油气藏的位置和大小。
全波形反演依赖于合成地震数据和实际地震数据的波形对比，目前在全波形反演中广泛采用的数据对比准则是最小二乘误差函数。然而当合成地震数据与实际地震数据的误差较大时，最小二乘误差函数易使全波形反演陷入局部极小值，从而导致最后得到的地下介质速度模型与真实模型相差较大，反演效果较差。
发明内容
本发明实施例提供一种全波形反演方法、设备及存储介质，用以解决现有基于最小二乘误差函数的全波形反演方法易陷入局部极小值，进而导致反演效果差的问题。
第一方面，本发明实施例提供一种全波形反演方法，包括：
步骤一：获取观测地震数据和初始速度模型；
步骤二：根据初始速度模型确定合成地震数据和正传波场；
步骤三：根据观测地震数据和合成地震数据构建匹配代价矩阵，匹配代价矩阵中的每一个矩阵元素用于表示合成地震数据中一个地震道数据与观测地震数据中一个地震道数据之间的相似性；
步骤四：根据匹配代价矩阵确定用于反映观测地震数据和合成地震数据误差的目标函数；
步骤五：计算目标函数关于合成地震数据的一阶偏导数，得到伴随场震源；
步骤六：根据伴随场震源与正传波场确定梯度；
步骤七：根据梯度和预设迭代步长对初始速度模型进行更新；
步骤八：重复执行步骤二至步骤七，直至满足预设收敛准则，将此时的速度模型作为全波形反演的最终速度模型输出。
一种实施例中，匹配代价矩阵中的矩阵元素根据如下表达式确定：
其中，cij表示匹配代价矩阵中第i行第j列的矩阵元素，cosθ表示pi和dj之间归一化的向量点积，pi表示合成地震数据在第i个检波器记录到的地震波形数据，dj表示观测地震数据在第j个检波器记录到的地震波形数据，λ≥0表示空间坐标权重，xi表示合成地震数据第i个检波器的空间位置坐标，xj表示观测地震数据第j个检波器的空间位置坐标，Ax表示归一化常数，用于对空间位置坐标进行归一化，Ai和Aj表示振幅权重，Ai根据观测地震数据第i个检波器记录到的地震波形确定，Aj根据观测地震数据第j个检波器记录到的地震波形确定。
一种实施例中，空间坐标权重与全波形反演各反演阶段的最大频率正相关。
一种实施例中，当用于声波全波形反演时，Ai＝Aj＝1。
一种实施例中，根据匹配代价矩阵确定用于反映观测地震数据和合成地震数据误差的目标函数包括：
采用竞价拍卖算法确定匹配代价矩阵的最优匹配路径；
将匹配代价矩阵中的位于最优匹配路径上的各矩阵元素的和确定为目标函数。
一种实施例中，目标函数的确定过程包括：
根据如下表达式确定使目标函数取值最小的矩阵M：
其中，Γ表示目标函数，mij表示M中第i行第j列的矩阵元素，i，j表示遍历匹配代价矩阵和矩阵M中的矩阵元素；
将矩阵M等价为置换向量σ，以使在匹配的过程中没有分裂，σi＝j表示将第i个元素映射到第j个元素，σi表示置换向量σ中第i个元素；
采用竞价拍卖算法进行求解，得到最优置换向量σ*；
根据最优置换向量σ*，采用如下表达式确定目标函数：
其中，表示最优置换向量σ*中第i个元素，/>表示匹配代价矩阵中第i行的元素，第/>列。
一种实施例中，伴随场震源根据如下表达式确定：
其中，si表示伴随场震源，表示向量pi和向量/>的内积，p表示合成地震数据中对应向量的模，d表示观测地震数据中对应向量的模。
一种实施例中，当观测地震数据包含的检波器数量小于合成地震数据包含的检波器数量时，对观测地震数据进行补零操作，使得观测地震数据包含的检波器数量等于合成地震数据包含的检波器数量。
第二方面，本发明实施例提供一种全波形反演设备，包括：
至少一个处理器和存储器；
存储器存储计算机执行指令；
至少一个处理器执行存储器存储的计算机执行指令，使得至少一个处理器执行如第一方面任一项所述的全波形反演方法。
第三方面，本发明实施例提供一种计算机可读存储介质，所述计算机可读存储介质中存储有计算机执行指令，计算机执行指令被处理器执行时用于实现如第一方面任一项所述的全波形反演方法。
本发明实施例提供的全波形反演方法、设备及存储介质，通过获取观测地震数据和初始速度模型，根据初始速度模型确定合成地震数据和正传波场，根据观测地震数据和合成地震数据构建匹配代价矩阵，匹配代价矩阵中的每一个矩阵元素用于表示合成地震数据中一个地震道数据与观测地震数据中一个地震道数据之间的相似性；根据匹配代价矩阵确定用于反映观测地震数据和合成地震数据误差的目标函数；计算目标函数关于合成地震数据的一阶偏导数，得到伴随场震源；根据伴随场震源与正传波场确定梯度；根据梯度和预设迭代步长对初始速度模型进行迭代更新，直至满足预设收敛准则，获得全波形反演的最终速度模型。将合成地震数据的地震道与观测地震数据的地震道进行匹配，充分考虑了地震道的空间特性，降低了对于初始速度模型的要求，不仅能够避免陷入局部极小值，而且能够提高反演精度。
附图说明
此处的附图被并入说明书中并构成本说明书的一部分，示出了符合本发明的实施例，并与说明书一起用于解释本发明的原理。
图1为本发明一实施例提供的全波形反演方法的流程图；
图2为本发明一实施例提供的透射模型的示意图；
图3为本发明一实施例提供的目标函数误差的示意图；
图4为本发明一实施例提供的盲测地震数据的示意图；
图5为采用现有全波形反演方法获得的反演结果的示意图；
图6为采用本发明实施例所提供的全波形反演方法获得的反演结果的示意图；
图7为本发明一实施例提供的全波形反演设备的结构示意图。
通过上述附图，已示出本发明明确的实施例，后文中将有更详细的描述。这些附图和文字描述并不是为了通过任何方式限制本发明构思的范围，而是通过参考特定实施例为本领域技术人员说明本发明的概念。
具体实施方式
下面通过具体实施方式结合附图对本发明作进一步详细说明。其中不同实施方式中类似元件采用了相关联的类似的元件标号。在以下的实施方式中，很多细节描述是为了使得本申请能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他元件、材料、方法所替代。在某些情况下，本申请相关的一些操作并没有在说明书中显示或者描述，这是为了避免本申请的核心部分被过多的描述所淹没，而对于本领域技术人员而言，详细描述这些相关操作并不是必要的，他们根据说明书中的描述以及本领域的一般技术知识即可完整了解相关操作。
另外，说明书中所描述的特点、操作或者特征可以以任意适当的方式结合形成各种实施方式。同时，方法描述中的各步骤或者动作也可以按照本领域技术人员所能显而易见的方式进行顺序调换或调整。因此，说明书和附图中的各种顺序只是为了清楚描述某一个实施例，并不意味着是必须的顺序，除非另有说明其中某个顺序是必须遵循的。
本文中为部件所编序号本身，例如“第一”、“第二”等，仅用于区分所描述的对象，不具有任何顺序或技术含义。而本申请所说“连接”、“联接”，如无特别说明，均包括直接和间接连接(联接)。
实施例一
图1为本发明一实施例提供的全波形反演方法的流程图。如图1所示，本实施例提供的全波形反演方法可以包括：
S101、获取观测地震数据和初始速度模型。
本实施例中的观测地震数据是实际采集的地震数据集，可以按照炮的编号进行存储。每一炮可以包含很多地震道，每一个地震道对应一个接收器。在野外的地震勘探中，炮的个数可以达到万的级别，每一炮也可以包含上万个地震道。本发明实施例所提供的全波形反演方法可以针对每一炮数据进行。本实施例中的初始速度模型可以采用现有方法获取，本实施例对此不作限制。
S102、根据初始速度模型确定合成地震数据和正传波场。
根据输入的初始速度模型，计算合成地震数据和正传波场。每一炮(炮也称为地震震源)都独立进行模拟。初始速度模型代表的是地下地质结构，具有一定的长度和深度。接收器(检波器)是布设在地表附近的数据记录仪。也就是说，对于每一炮，都会有其对应的正传波场以及合成地震数据集合。正传波场是在整个初始速度模型空间的波场，合成地震数据集合指的是对应于接收器位置的波场记录。通常在地震行业，无论是二维模型还是三维模型，计算的合成地震数据都可以存储在一个二维矩阵里面，每一列对应一个检波器记录到的地震道信号，每一行都是按照接收器的编号存储。
S103、根据观测地震数据和合成地震数据构建匹配代价矩阵，匹配代价矩阵中的每一个矩阵元素用于表示合成地震数据中一个地震道数据与观测地震数据中一个地震道数据之间的相似性。
本实施例中可以针对每一个炮集构建匹配代价矩阵，这是因为对于一个给定的炮号，则与之相对应的有一个观测地震数据炮集和一个合成地震数据炮集。匹配代价矩阵中的每一个矩阵元素对应一对地震道，一个合成数据道集，一个观测数据道集。
S104、根据匹配代价矩阵确定用于反映观测地震数据和合成地震数据误差的目标函数。
可以在匹配代价矩阵中，寻找一个最优化路径，使得沿着该路径的矩阵元素的求和最优。该最优化路径反映了观测地震数据和合成地震数据中相匹配的地震道集，用匹配的地震道集计算目标函数。
S105、计算目标函数关于合成地震数据的一阶偏导数，得到伴随场震源。
S106、根据伴随场震源与正传波场确定梯度。
本实施例中在得到伴随场震源之后，可以反传伴随场震源。将伴随场与正传波场相乘，并累加求和，得到梯度。需要说明的是，此处的正传波场通常有两种获得方式：一是在步骤S102中做正传波场计算的时候，把整个正传波场存储在硬盘上，然后在步骤S106中从硬盘中把正传波场读进来；二是在步骤S102中进行正传波场计算的时候，记录并存储边界值，而并不需要将整个正传波场存储在硬盘上，然后在步骤S106中根据存储的边界值重构整个地震正传波场。
S107、根据梯度和预设迭代步长对初始速度模型进行更新。重复执行步骤S102至步骤S107，直至满足预设收敛准则，将此时的速度模型作为全波形反演的最终速度模型输出。
在上一步求出梯度后，这一步根据梯度对初始速度模型进行迭代更新。本实施例中的预设迭代步长可以是事先给定的一个固定迭代步长，也可以是在迭代更新过程中根据目标函数的取值不断进行优化的优化迭代步长。
在对初始速度模型进行更新之后，会获得新模型。使用新模型替换原有初始速度模型，然后重复执行步骤S102至步骤S107，直至满足预设收敛准则，得到最终速度模型。预设收敛准则例如可以是观测地震数据与合成地震数据之间的误差小于预设阈值。
本实施例提供的全波形反演方法，通过获取观测地震数据和初始速度模型，根据初始速度模型确定合成地震数据和正传波场，根据观测地震数据和合成地震数据构建匹配代价矩阵，匹配代价矩阵中的每一个矩阵元素用于表示合成地震数据中一个地震道数据与观测地震数据中一个地震道数据之间的相似性；根据匹配代价矩阵确定用于反映观测地震数据和合成地震数据误差的目标函数；计算目标函数关于合成地震数据的一阶偏导数，得到伴随场震源；根据伴随场震源与正传波场确定梯度；根据梯度和预设迭代步长对初始速度模型进行迭代更新，直至满足预设收敛准则，获得全波形反演的最终速度模型。将合成地震数据的地震道与观测地震数据的地震道进行匹配，充分考虑了地震道的空间特性，降低了对于初始速度模型的要求，不仅能够避免陷入局部极小值，而且能够提高反演精度。
实施例二
在上述实施例的基础上，本实施例进一步地详细阐述如何构建匹配代价矩阵。在地震勘探里，检波器通常可以看作一个点。在匹配合成地震数据集和实际观测地震数据集的时候，可以把检波器的位置看作是锚点，把地震信号和位置信息看作是锚点的属性向量。可以采用dj表示观测地震数据在第j个检波器记录到的地震波形数据；采用符号pi表示合成地震数据在第i个检波器记录到的地震波形数据。若将观测地震数据包含的检波器数量记为n，则j≤n；若将合成地震数据包含的检波器数量记为m，则i≤m。在二维平地表情形下，每一个检波器用一个坐标就可以表征。首先建立一个匹配代价矩阵，这个矩阵的每一个元素代表两个地震道之间的相似性。
一种可选的实施方式中，匹配代价矩阵中的矩阵元素可以根据如下表达式确定：
其中，cij表示匹配代价矩阵中第i行第j列的矩阵元素，cosθ表示pi和dj之间归一化的向量点积，θ表示pi和dj之间的向量夹角，pi表示合成地震数据在第i个检波器记录到的地震波形数据，dj表示观测地震数据在第j个检波器记录到的地震波形数据，λ≥0表示空间坐标权重，xi表示合成地震数据第i个检波器的空间位置坐标，xj表示观测地震数据第j个检波器的空间位置坐标，Ax表示归一化常数，用于对空间位置坐标进行归一化，可以取最大的空间位置坐标数值。Ai和Aj表示振幅权重，Ai根据观测地震数据第i个检波器记录到的地震波形确定，Aj根据观测地震数据第j个检波器记录到的地震波形确定。一种可选的实施方式中，当用于声波全波形反演时，可以选用Ai＝Aj＝1。可选的，也可以采用均方根来创建振幅函数。由于采用合成地震数据的波形来创建振幅权重会导致伴随场震源变得比较复杂，为了避免出现这样的问题，本实施例中采用观测地震数据的波形来创建振幅权重。
一种可选的实施方式中，空间坐标权重与全波形反演各反演阶段的最大频率正相关。空间坐标权重决定了最优化的置换向量。考虑以下两种极端情况：一种极端情况下，空间坐标权重的取值是λ＝0，这时候空间位置坐标变量完全不起作用，最优化匹配只取决于波形的相似度；另一种极端情形下λ＝∞，波形的振幅可以忽略，最优化匹配函数等价于归一化的最小二乘或者零延迟的互相关函数。选择一个较小的空间坐标权重λ，那么全波形反演的目标函数会更凸，可以适当解决陷入局部极小值的问题；而一个较大的空间坐标权重λ会使全波形反演的目标函数更接近于最小二乘，这对于提高反演精度很有用。因此在全波形反演中，可以在后面的反演阶段逐渐增大空间坐标权重λ。也就是说，空间坐标权重λ可以随着反演阶段的最大频率的增大而增大，即空间坐标权重与全波形反演各反演阶段的最大频率正相关，这样既可以避免陷入局部极小值，同时又可以重构出精细的地下地质结构。
一种可选的实施例中，根据匹配代价矩阵确定用于反映观测地震数据和合成地震数据误差的目标函数包括：采用竞价拍卖算法确定匹配代价矩阵的最优匹配路径；将匹配代价矩阵中的位于最优匹配路径上的各矩阵元素的和确定为目标函数。
具体的，目标函数的确定过程如下：
首先根据如下表达式确定使目标函数取值最小的矩阵M：
其中，Γ表示目标函数，mij表示M中第i行第j列的矩阵元素，i，j表示遍历匹配代价矩阵和矩阵M中的矩阵元素。mij只能取0或者1，同时为了保证每一个点在匹配过程中没有分裂，还需要限制矩阵M每一行以及每一列都只能有一个1。
可以将其看作是一个线性分配问题，将矩阵M等价为置换向量σ，以使在匹配的过程中没有分裂，σi＝j表示将第i个元素映射到第j个元素，或者说将元素i映射到元素j，σi表示置换向量σ中第i个元素。
采用竞价拍卖算法进行求解，得到最优置换向量σ*。通过竞价拍卖算法可以快速求解得到一个最优化的置换向量σ*。
那么总的匹配代价，也就是全波形反演里要用到的目标函数，可以根据最优置换向量σ*，采用如下表达式来确定：
其中，表示最优置换向量σ*中第i个元素，/>表示匹配代价矩阵中第i行，第/>列的元素。
在求解出最优化的置换向量和目标函数值之后，下一步便是求解伴随场震源。伴随场震源si是目标函数Γ关于合成地震数据pi的一阶偏导数，一种可选的实施方式中，伴随场震源可以根据如下表达式确定：
其中，si表示伴随场震源，表示向量pi和向量/>的内积，p表示合成地震数据中对应向量的模(二范数下的长度)，d表示观测地震数据中对应向量的模(二范数下的长度)。/>表示观测地震数据中与合成地震数据中第i个检波器记录到的地震波形数据相匹配的地震波形数据。
本发明实施例所提供的全波形反演方法对于非均匀分布的检波器也是适用的，因为匹配代价矩阵把空间变量也包括进去了。本发明实施例所提供的全波形反演方法对于三维数据采集也是适用的，针对三维数据的情形，只需要在构建匹配代价矩阵的时候把另一维度的空间变量也加进去即可。具体地，可以先把所有的检波器形成的锚点用一维数组存储，然后在构建匹配代价矩阵的时候把两个空间变量和振幅都包含进去，这样形成的还是一个二维的匹配代价矩阵，后续处理过程依然可以使用同一套方法。
实施例三
本实施例中采用具体测试数据来对本发明实施例所提供的全波形反演方法进行进一步说明。请参考图2，图2为本发明一实施例提供的透射模型的示意图。如图2所示的透射模型，400个接收器在下边界，单个震源在上边界中间。震源子波是5赫兹的Ricker子波。左右两个速度异常都是根据高斯函数生成，左侧为低速异常，右侧为高速异常。最低速和最高速是两个独立的参数，共计生成51X51＝2601个速度模型。然后对于每一个速度模型，计算地震记录。基于图2所示的透射模型下的目标函数误差图如图3所示，图3为本发明一实施例提供的目标函数误差的示意图。选取一个速度模型，然后计算所有其他的地震记录跟参考记录的误差。图3中左侧示意图是用传统方法，也就是最小二乘误差函数得到的，可以看出可行区域很狭窄，对初始模型要求较高。图3中右侧示意图是基于本发明实施例所提供的最优化匹配代价矩阵得到的，凸性良好。图4为本发明一实施例提供的盲测地震数据的示意图。图4为Chevron 2014的盲测地震数据的示意图，图4中所展示的是20赫兹以下的数据。噪声在低频段很高，对于全波形反演来说挑战很大。图5为采用现有全波形反演方法获得的反演结果的示意图。如图5所示，当采用常规的全波形反演方法时，反演的模型出现速度异常，中间2到3公里深度的低速层不清晰，这是由于地震数据的反射波没有匹配好。图6为采用本发明实施例所提供的全波形反演方法获得的反演结果的示意图。如图6所示，采用本发明实施例所提供的全波形反演方法获得的反演结果，反射层位清晰，速度结构合理。进一步验证了本发明实施例所提供方法的能够提高反演效果。
实施例四
通常情况下合成地震数据所包含的检波器数量和位置跟观测地震数据一一对应，从而合成地震数据与实际观测地震数据具有一样的尺寸，即两者具有同样的地震道个数，同样的记录时间，采用本发明实施例所提供的全波形反演方法可以快速准确的获得高精度的反演结果，既可以实现平衡性匹配。
在全波形反演的各种实际应用场景中，有时候检波器并不是密集分布，或者说分布较稀疏且不规则，也就是说实际应用场景下获得的观测地震数据可能较为稀疏。为了实现在一个覆盖更密集的观测系统里求解全波形反演问题，可以在合成地震数据中设置更多的检波器。此时合成地震数据比观测地震数据更加密集。若采用平衡性匹配算法的话，检波器将会被匹配到一个较远的距离。为了解决这个问题，可以采用冗余性匹配算法。实际地震数据的炮集包含n个检波器，在全波形反演中，设置m个接收器，通常m>n以保证覆盖较为密集。此时需要将n个点匹配到m个点上，可以把观测地震数据补齐到m个检波器，把这些新添加的检波器的地震信号设为零。最终，将冗余性匹配转换为平衡性匹配，可以采用平衡性匹配算法来实现，这样做的好处是匹配算法可以在一个覆盖更密集的观测系统里求解全波形反演问题。
也就是说，在一种可选的实施方式中，当观测地震数据包含的检波器数量小于合成地震数据包含的检波器数量时，对观测地震数据进行补零操作，使得观测地震数据包含的检波器数量等于合成地震数据包含的检波器数量。
实施例五
本发明实施例还提供一种全波形反演设备，请参见图7所示，本发明实施例仅以图7为例进行说明，并不表示本发明仅限于此。图7为本发明一实施例提供的全波形反演设备的结构示意图。如图7所示，本实施例提供的全波形反演设备70包括：存储器701、处理器702和总线703。其中，总线703用于实现各元件之间的连接。
存储器701中存储有计算机程序，计算机程序被处理器702执行时可以实现上述任一方法实施例的技术方案。
其中，存储器701和处理器702之间直接或间接地电性连接，以实现数据的传输或交互。例如，这些元件相互之间可以通过一条或者多条通信总线或信号线实现电性连接，如可以通过总线703连接。存储器701中存储有实现全波形反演方法的计算机程序，包括至少一个可以软件或固件的形式存储于存储器701中的软件功能模块，处理器702通过运行存储在存储器701内的软件程序以及模块，从而执行各种功能应用以及数据处理。
存储器701可以是，但不限于，随机存取存储器(Random Access Memory，简称：RAM)，只读存储器(Read Only Memory，简称：ROM)，可编程只读存储器(ProgrammableRead-Only Memory，简称：PROM)，可擦除只读存储器(Erasable Programmable Read-OnlyMemory，简称：EPROM)，电可擦除只读存储器(Electric Erasable Programmable Read-Only Memory，简称：EEPROM)等。其中，存储器701用于存储程序，处理器702在接收到执行指令后，执行程序。进一步地，上述存储器701内的软件程序以及模块还可包括操作系统，其可包括各种用于管理系统任务(例如内存管理、存储设备控制、电源管理等)的软件组件和/或驱动，并可与各种硬件或软件组件相互通信，从而提供其他软件组件的运行环境。
处理器702可以是一种集成电路芯片，具有信号的处理能力。上述的处理器702可以是通用处理器，包括中央处理器(Central Processing Unit，简称：CPU)、网络处理器(Network Processor，简称：NP)等。可以实现或者执行本发明实施例中的公开的各方法、步骤及逻辑框图。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。可以理解，图7的结构仅为示意，还可以包括比图7中所示更多或者更少的组件，或者具有与图7所示不同的配置。图7中所示的各组件可以采用硬件和/或软件实现。
本发明实施例还提供一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行以实现上述任一方法实施例的技术方案。
本公开中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。
本公开的保护范围不限于上述的实施例，显然，本领域的技术人员可以对本公开进行各种改动和变形而不脱离本公开的范围和精神。倘若这些改动和变形属于本公开权利要求及其等同技术的范围，则本公开的意图也包含这些改动和变形在内。
Claims (10)
1.一种全波形反演方法，其特征在于，包括：
步骤一：获取观测地震数据和初始速度模型；
步骤二：根据所述初始速度模型确定合成地震数据和正传波场；
步骤三：根据所述观测地震数据和所述合成地震数据构建匹配代价矩阵，所述匹配代价矩阵中的每一个矩阵元素用于表示所述合成地震数据中一个地震道数据与所述观测地震数据中一个地震道数据之间的相似性；
步骤四：根据所述匹配代价矩阵确定用于反映所述观测地震数据和所述合成地震数据误差的目标函数；
步骤五：计算所述目标函数关于所述合成地震数据的一阶偏导数，得到伴随场震源；
步骤六：根据所述伴随场震源与所述正传波场确定梯度；
步骤七：根据所述梯度和预设迭代步长对所述初始速度模型进行更新；
步骤八：重复执行步骤二至步骤七，直至满足预设收敛准则，将此时的速度模型作为全波形反演的最终速度模型输出。
2.根据权利要求1所述的方法，其特征在于，所述匹配代价矩阵中的矩阵元素根据如下表达式确定：
其中，cij表示匹配代价矩阵中第i行第j列的矩阵元素，cosθ表示pi和dj之间归一化的向量点积，pi表示合成地震数据在第i个检波器记录到的地震波形数据，dj表示观测地震数据在第j个检波器记录到的地震波形数据，λ≥0表示空间坐标权重，xi表示合成地震数据第i个检波器的空间位置坐标，xj表示观测地震数据第j个检波器的空间位置坐标，Ax表示归一化常数，用于对空间位置坐标进行归一化，Ai和Aj表示振幅权重，Ai根据观测地震数据第i个检波器记录到的地震波形确定，Aj根据观测地震数据第j个检波器记录到的地震波形确定。
3.根据权利要求2所述的方法，其特征在于，所述空间坐标权重与全波形反演各反演阶段的最大频率正相关。
4.根据权利要求2所述的方法，其特征在于，当用于声波全波形反演时，Ai＝Aj＝1。
5.根据权利要求2所述的方法，其特征在于，所述根据所述匹配代价矩阵确定用于反映所述观测地震数据和所述合成地震数据误差的目标函数包括：
采用竞价拍卖算法确定所述匹配代价矩阵的最优匹配路径；
将所述匹配代价矩阵中的位于所述最优匹配路径上的各矩阵元素的和确定为所述目标函数。
6.根据权利要求5所述的方法，其特征在于，所述目标函数的确定过程包括：
根据如下表达式确定使目标函数取值最小的矩阵M：
其中，Γ表示目标函数，mij表示M中第i行第j列的矩阵元素，i,j表示遍历匹配代价矩阵和矩阵M中的矩阵元素；
将矩阵M等价为置换向量σ，以使在匹配的过程中没有分裂，σi＝j表示将第i个元素映射到第j个元素，σi表示置换向量σ中第i个元素；
采用竞价拍卖算法进行求解，得到最优置换向量σ*；
根据所述最优置换向量σ*，采用如下表达式确定所述目标函数：
其中，表示所述最优置换向量σ*中第i个元素，/>表示所述匹配代价矩阵中第i行，第/>列的元素。
7.根据权利要求6所述的方法，其特征在于，所述伴随场震源根据如下表达式确定：
其中，si表示伴随场震源，表示向量pi和向量/>的内积，p表示合成地震数据中对应向量的模，d表示观测地震数据中对应向量的模。
8.根据权利要求1-7任一项所述的方法，其特征在于，当所述观测地震数据包含的检波器数量小于所述合成地震数据包含的检波器数量时，对所述观测地震数据进行补零操作，使得所述观测地震数据包含的检波器数量等于所述合成地震数据包含的检波器数量。
9.一种全波形反演设备，其特征在于，包括：至少一个处理器和存储器；
所述存储器存储计算机执行指令；
所述至少一个处理器执行所述存储器存储的计算机执行指令，使得所述至少一个处理器执行如权利要求1-8任一项所述的全波形反演方法。
10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质中存储有计算机执行指令，所述计算机执行指令被处理器执行时用于实现如权利要求1-8任一项所述的全波形反演方法。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211105068.7A CN117687083A (zh) | 2022-09-09 | 2022-09-09 | 全波形反演方法、设备及存储介质 |
| PCT/CN2023/117802 WO2024051834A1 (zh) | 2022-09-09 | 2023-09-08 | 全波形反演方法、设备及存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211105068.7A CN117687083A (zh) | 2022-09-09 | 2022-09-09 | 全波形反演方法、设备及存储介质 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117687083A true CN117687083A (zh) | 2024-03-12 |
Family
ID=90132567
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211105068.7A Pending CN117687083A (zh) | 2022-09-09 | 2022-09-09 | 全波形反演方法、设备及存储介质 |
Country Status (2)
| Country | Link |
|---|---|
| CN (1) | CN117687083A (zh) |
| WO (1) | WO2024051834A1 (zh) |
Family Cites Families (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103472482B (zh) * | 2013-09-03 | 2016-07-06 | 中国石油天然气集团公司 | 基于基因排序体系的多波地震资料时间域匹配方法及系统 |
| WO2016001750A2 (en) * | 2014-06-30 | 2016-01-07 | Cgg Services Sa | Seismic data processing using matching filter based cost function optimization |
| GB2538807B (en) * | 2015-05-29 | 2019-05-15 | Sub Salt Solutions Ltd | Method for improved geophysical investigation |
| GB2538804A (en) * | 2015-05-29 | 2016-11-30 | Sub Salt Solutions Ltd | Improved method for inversion modelling |
| CN111045077B (zh) * | 2019-12-20 | 2022-03-18 | 核工业北京地质研究院 | 一种陆地地震数据的全波形反演方法 |
| CN113970789B (zh) * | 2020-07-24 | 2024-04-09 | 中国石油化工股份有限公司 | 全波形反演方法、装置、存储介质及电子设备 |
| CN114966837A (zh) * | 2022-05-13 | 2022-08-30 | 中国石油大学(北京) | 基于卷积型w2距离目标函数的全波形反演方法及装置 |
-
2022
- 2022-09-09 CN CN202211105068.7A patent/CN117687083A/zh active Pending
-
2023
- 2023-09-08 WO PCT/CN2023/117802 patent/WO2024051834A1/zh unknown
Also Published As
| Publication number | Publication date |
|---|---|
| WO2024051834A1 (zh) | 2024-03-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| AU2016331881B8 (en) | Q-compensated full wavefield inversion | |
| KR101797451B1 (ko) | 상호상관 목적 함수를 통한 해양 스트리머 데이터에 대한 동시 소스 반전 | |
| EP3063562B1 (en) | Methods of subsurface exploration, computer program product and computer-readable storage medium | |
| US9632192B2 (en) | Method of processing seismic data by providing surface offset common image gathers | |
| US10557956B2 (en) | Method and system of processing seismic data by providing surface aperture common image gathers | |
| WO2004090573A2 (en) | Seimsic imaging by wave migration using a krylov space expansion of the square root exponent operator | |
| WO2009002001A1 (en) | Method for velocity analysis using waveform inversion in laplace domain for geophysical imaging | |
| US10310117B2 (en) | Efficient seismic attribute gather generation with data synthesis and expectation method | |
| US20160161620A1 (en) | System and method for estimating repeatability using base data | |
| WO2013093468A2 (en) | Full waveform inversion quality control method | |
| RU2339056C2 (ru) | Обобщенное трехмерное прогнозирование кратных волн от поверхности | |
| WO2016193180A1 (en) | Improved method for inversion modelling | |
| WO2013093467A1 (en) | Method of, and apparatus for, full waveform inversion | |
| BR102013004904A2 (pt) | Métodos e aparelho para remoção automatizada de ruído de dados sísmicos | |
| Thiel et al. | Comparison of acoustic and elastic full‐waveform inversion of 2D towed‐streamer data in the presence of salt | |
| US9921324B2 (en) | Systems and methods employing upward beam propagation for target-oriented seismic imaging | |
| GB2481270A (en) | Method of, and apparatus for, full waveform inversion modelling | |
| US11397273B2 (en) | Full waveform inversion in the midpoint-offset domain | |
| EP3092513A2 (en) | Systems and methods for destriping seismic data | |
| GB2503640A (en) | Quality Assurance in a Full Waveform Inversion Process | |
| CN117687083A (zh) | 全波形反演方法、设备及存储介质 | |
| US20240159930A1 (en) | Method and apparatus for implementing full waveform inversion using angle gathers | |
| US20230051004A1 (en) | Method and apparatus for performing efficient modeling of extended-duration moving seismic sources | |
| WO2016071728A1 (en) | Systems and methods for vortex calculation as attribute for geologic discontinuities |
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 |