CN116629031A - 一种全波形反演方法、装置、电子设备及存储介质 - Google Patents
一种全波形反演方法、装置、电子设备及存储介质 Download PDFInfo
- Publication number
- CN116629031A CN116629031A CN202310893372.0A CN202310893372A CN116629031A CN 116629031 A CN116629031 A CN 116629031A CN 202310893372 A CN202310893372 A CN 202310893372A CN 116629031 A CN116629031 A CN 116629031A
- Authority
- CN
- China
- Prior art keywords
- operator
- compensation
- source
- determining
- seismic source
- 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 75
- 238000004088 simulation Methods 0.000 claims abstract description 35
- 230000036962 time dependent Effects 0.000 claims abstract description 20
- 230000006870 function Effects 0.000 claims description 61
- 239000006185 dispersion Substances 0.000 claims description 20
- 230000005540 biological transmission Effects 0.000 claims description 11
- 238000013016 damping Methods 0.000 claims description 6
- 230000001419 dependent effect Effects 0.000 claims description 6
- 238000010586 diagram Methods 0.000 description 15
- 230000008569 process Effects 0.000 description 14
- 238000004891 communication Methods 0.000 description 8
- 230000006641 stabilisation Effects 0.000 description 8
- 238000011105 stabilization Methods 0.000 description 8
- 238000004364 calculation method Methods 0.000 description 7
- 238000004590 computer program Methods 0.000 description 7
- 230000009471 action Effects 0.000 description 5
- 230000008901 benefit Effects 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 230000003287 optical effect Effects 0.000 description 3
- 230000036961 partial effect Effects 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 2
- 238000005286 illumination Methods 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 230000000704 physical effect Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 238000013473 artificial intelligence Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000001413 cellular effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 239000013307 optical fiber Substances 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 230000001953 sensory effect Effects 0.000 description 1
- 230000005477 standard model Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V20/00—Geomodelling in general
-
- 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/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- 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
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Operations Research (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本申请提供了一种全波形反演方法、装置、电子设备及存储介质；所述方法包括：设定震源，并构建所述震源对应的目标函数；基于所述震源对应的模拟数据，确定所述震源对应的震源波场；基于所述震源对应的模拟数据和观测数据，确定所述震源对应的逆时伴随波场；基于所述震源波场和所述逆时伴随波场，确定所述目标函数的补偿梯度；根据所述补偿梯度对所述目标函数的初始参数进行迭代更新，直至所述迭代更新的次数达到预设值。如此，能够智能地进行全波形反演，提高了全波形反演的精度和效率。
Description
技术领域
本申请涉及全波形反演技术，尤其涉及一种全波形反演方法、装置、电子设备及存储介质。
背景技术
现有的大多数全波形反演是使用流变模型作为控制方程，由于流变模型表征的模型参数相对较多，且振幅衰减和相位频散作用相互耦合，现有的全波形反演的精度较低，且全波形反演的计算过程的十分耗时。人们更希望减少全波形反演的计算过程的时间并提高全波形反演的精度。
因此，如何智能地进行全波形反演，以提高全波形反演的精度和效率是一直追求的目标。
发明内容
本申请实施例提供了一种全波形反演方法、装置、电子设备及存储介质。
根据本申请的第一方面，提供了一种全波形反演方法，该方法包括：设定震源，并构建所述震源对应的目标函数；基于所述震源对应的模拟数据，确定所述震源对应的震源波场；基于所述震源对应的模拟数据和观测数据，确定所述震源对应的逆时伴随波场；基于所述震源波场和所述逆时伴随波场，确定所述目标函数的补偿梯度；根据所述补偿梯度对所述目标函数的初始参数进行迭代更新，直至所述迭代更新的次数达到预设值。
根据本申请一实施方式，所述基于所述震源对应的模拟数据，确定所述震源对应的震源波场，包括：基于所述模拟数据，确定对应的第一稳定补偿算子；所述第一稳定补偿算子包括振幅补偿算子、高频稳定算子、相位频散算子、相位校正算子中的至少之一；基于所述第一稳定补偿算子，确定所述模拟数据对应的波场延拓算子；基于所述波场延拓算子和所述震源，以时间正向的方式求解补偿方程，得的所述震源对应的正传波场；对所述正传波场以时间反向的方式进行重构，得到震源波场。
根据本申请一实施方式，所述基于所述震源对应的模拟数据和观测数据，确定所述震源对应的逆时伴随波场，包括：基于所述模拟数据，确定对应的第二稳定补偿算子；所述第二稳定补偿算子包括振幅补偿算子、相位频散算子中的至少之一；基于所述第二稳定补偿算子，求解衰减方程，得到所述震源对应的粘声波场；基于所述粘声波场和所述观测数据，确定伴随震源；基于所述伴随震源，确定所述逆时伴随波场。
根据本申请一实施方式，所述基于所述伴随震源，确定所述逆时伴随波场，包括：基于所述伴随震源，确定对应的第三稳定补偿算子；所述第三稳定补偿算子包括振幅补偿算子、高频稳定算子、相位频散算子、相位校正算子中的至少之一；基于所述第三稳定补偿算子，确定所述伴随震源对应的伴随算子；基于所述伴随算子和所述伴随震源，以时间反向的方式求解伴随方程，得到所述震源对应的逆时伴随波场。
根据本申请一实施方式，所述基于所述震源波场和所述逆时伴随波场，确定所述目标函数的补偿梯度，包括：将所述震源波场和所述逆时伴随波场互相关，得到所述目标函数的补偿梯度。
根据本申请一实施方式，所述根据所述补偿梯度对所述目标函数的初始参数进行迭代更新，直至所述迭代更新的次数达到预设值之后，所述方法还包括：将最后一次迭代更新后的参数作为所述目标函数对应的全波形反演结果。
根据本申请一实施方式，所述振幅补偿算子、高频稳定算子、相位频散算子、相位校正算子均包括固定阶数的拉普拉斯算子。
根据本申请的第二方面，提供了一种全波形反演装置，所述全波形反演装置包括：设定模块，用于设定震源，并构建所述震源对应的目标函数；第一确定模块，用于基于所述震源对应的模拟数据，确定所述震源对应的震源波场；第二确定模块，用于基于所述震源对应的模拟数据和观测数据，确定所述震源对应的逆时伴随波场；梯度确定模块，用于基于所述震源波场和所述逆时伴随波场，确定所述目标函数的补偿梯度；迭代模块，用于根据所述补偿梯度对所述目标函数的初始参数进行迭代更新，直至所述迭代更新的次数达到预设值。
根据本申请的第三方面，提供了一种电子设备，包括：
至少一个处理器；以及
与所述至少一个处理器通信连接的存储器；其中，
所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被所述至少一个处理器执行，以使所述至少一个处理器能够执行本申请所述的方法。
根据本申请的第四方面，提供了一种存储有计算机指令的非瞬时计算机可读存储介质，所述计算机指令用于使所述计算机执行本申请所述的方法。
本申请实施例的方法，设定震源，并构建所述震源对应的目标函数；基于所述震源对应的模拟数据，确定所述震源对应的震源波场；基于所述震源对应的模拟数据和观测数据，确定所述震源对应的逆时伴随波场；基于所述震源波场和所述逆时伴随波场，确定所述目标函数的补偿梯度；根据所述补偿梯度对所述目标函数的初始参数进行迭代更新，直至所述迭代更新的次数达到预设值。如此，能够智能地进行全波形反演，提高了全波形反演的精度和效率。
需要理解的是，本申请的教导并不需要实现上面所述的全部有益效果，而是特定的技术方案可以实现特定的技术效果，并且本申请的其他实施方式还能够实现上面未提到的有益效果。
附图说明
通过参考附图阅读下文的详细描述，本申请示例性实施方式的上述以及其他目的、特征和优点将变得易于理解。在附图中，以示例性而非限制性的方式示出了本申请的若干实施方式，其中：
在附图中，相同或对应的标号表示相同或对应的部分。
图1示出了本申请实施例提供的全波形反演方法的处理流程示意图一；
图2示出了本申请实施例提供的全波形反演方法的处理流程示意图二；
图3示出了本申请实施例提供的全波形反演方法的处理流程示意图三；
图4示出了本申请实施例提供的全波形反演方法的处理流程示意图四；
图5示出了本申请实施例提供的全波形反演方法的处理流程示意图五；
图6示出了本申请实施例提供的全波形反演方法的一种应用场景图；
图7示出了本申请实施例提供的全波形反演方法的另一种应用场景图；
图8示出了本申请实施例提供的全波形反演装置的一种可选示意图；
图9示出了本申请实施例提供的电子设备的组成结构示意图。
具体实施方式
为使本申请的目的、特征、优点能够更加的明显和易懂，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而非全部实施例。基于本申请中的实施例，本领域技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。
在以下的描述中，涉及到“一些实施例”，其描述了所有可能实施例的子集，但是可以理解，“一些实施例”可以是所有可能实施例的相同子集或不同子集，并且可以在不冲突的情况下相互结合。
对本申请实施例提供的全波形反演方法中的处理流程进行说明。参见图1，图1是本申请实施例提供的全波形反演方法的处理流程示意图一，将结合图1示出的步骤S101-S105进行说明。
步骤S101，设定震源，并构建震源对应的目标函数。
在一些实施例中，震源可以包括：雷克子波。震源还可以包括其他地震波，本申请实施例不作限定。目标函数可以包括：L2目标函数。本申请实施例不限定具体的目标函数。构建震源对应的目标函数可以包括：设定目标函数的初始参数，初始参数可以为速度参数。
震源对应的目标函数可以通过下述公式（1）表示：
（1）
其中，表示目标函数，m表示在空间x上的模型参数，d表示观测数据，p表示模拟数据，/>表示检波器位置，t表示时间。
步骤S102，基于震源对应的模拟数据，确定震源对应的震源波场。
在一些实施例中，步骤S102可以包括：基于模拟数据，确定对应的第一稳定补偿算子；基于第一稳定补偿算子，确定模拟数据对应的波场延拓算子；基于波场延拓算子和震源，以时间正向的方式求解补偿方程，得的震源对应的正传波场；对正传波场以时间反向的方式进行重构，得到震源波场。其中，第一稳定补偿算子可以包括振幅补偿算子、高频稳定算子、相位频散算子、相位校正算子中的至少之一；观测数据可以包括：采用多个检波点组成观测系统，通过观测系统获得观测数据。模拟数据可以包括建立速度场模型，并根据震源情况计算得出模拟数据。振幅补偿算子、高频稳定算子、相位频散算子、相位校正算子中均包括固定阶数的拉普拉斯算子。对正传波场以时间反向的方式进行重构，得到震源波场可以包括：将正传波场存储至存储器中，再从存储器中读取正传波场，得到震源波场。
针对基于模拟数据，确定对应的第一稳定补偿算子可以通过下述公式（2）~（7）表示：
（2）
（3）
（4）
（5）
（6）
（7）
其中，表示在参考频率/>下的相速度，c表示介质速度。/>表示模拟频率。/>表示拉普拉斯算子，Q表示品质因子。/>和/>表示稳定补偿参数，其数值通常取决于介质的物理性质，/>，/>。一般来说，较大的/>可以保留更多的高频信息，且随着/>的增大更多的高波数信息被压制。/>表示相位频散算子，/>表示振幅补偿算子，/>和/>表示相位校正算子，/>表示高频稳定算子。p表示模拟数据，t表示时间。
针对基于第一稳定补偿算子，确定模拟数据对应的波场延拓算子A可以通过下述公式（8）表示：
（8）
其中，表示相位频散算子，/>表示振幅补偿算子，/>和/>表示相位校正算子，/>表示高频稳定算子。/>，/>，/>表示模拟频率，/>表示参考频率，Q表示品质因子，p表示模拟数据，t表示时间。波场延拓算子A为描述地震波在地下传播的偏微分方程组（Partial Differential Equations，PDEs）。因此，全波形反演可以被看作是一个PDEs的约束优化问题。此外，该偏微分方程组具有一个独特的优势，即波场延拓算子A可解耦表示为不同的独立算子，包括振幅补偿算子、高频稳定算子、相位频散算子和相位校正算子。在实际计算中，可根据具体的需求组成运动学和动力学信息解耦的全波形反演梯度。
针对基于波场延拓算子A和震源f，以时间正向的方式求解补偿方程，得的震源对应的正传波场，可以通过下述公式（9）表示：
（9）
其中，表示正传波场，公式（9）为补偿方程。
步骤S103，基于震源对应的模拟数据和观测数据，确定震源对应的逆时伴随波场。
在一些实施例中，步骤S103可以包括：基于模拟数据，确定对应的第二稳定补偿算子；基于第二稳定补偿算子，求解衰减方程，得到震源对应的粘声波场；基于粘声波场和观测数据，确定伴随震源；基于伴随震源，确定逆时伴随波场。
针对基于模拟数据，确定对应的第二稳定补偿算子，可以通过上述公式（2）和（3）表示，第二稳定补偿算子可以包括振幅补偿算子、相位频散算子中的至少之一。
针对基于第二稳定补偿算子，求解衰减方程，得到震源对应的粘声波场，可以通过下述公式（10）表示：
（10）
其中，表示相位频散算子，/>表示振幅补偿算子，c表示介质速度，/>表示粘声波场，/>，/>表示模拟频率，/>表示参考频率，/>，Q表示品质因子。t表示时间。公式（10）为衰减方程。
针对基于粘声波场和观测数据，确定伴随震源，可以包括：将粘声波场与观测数据d相减，得到粘声波场与观测数据的残差，将该残差作为伴随震源。
针对基于伴随震源，确定逆时伴随波场，在具体实施时可以包括：基于伴随震源，确定对应的第三稳定补偿算子；基于第三稳定补偿算子，确定伴随震源对应的伴随算子；基于伴随算子和伴随震源，以时间反向的方式求解伴随方程，得到震源对应的逆时伴随波场。
针对基于伴随震源，确定对应的第三稳定补偿算子，可以通过下述公式（11）~（15）表示：
（11）
（12）
（13）
（14）
（15）
其中，c表示介质速度。表示模拟频率。/>表示拉普拉斯算子，Q表示品质因子。和/>表示稳定补偿参数，其数值通常取决于介质的物理性质，/>，/>。一般来说，较大的/>可以保留更多的高频信息，且随着/>的增大更多的高波数信息被压制。/>表示相位频散算子，/>表示振幅补偿算子，/>和/>表示相位校正算子，/>表示高频稳定算子。μ表示伴随波场，t表示时间。
针对基于第三稳定补偿算子，确定伴随震源对应的伴随算子，可以通过下述公式（16）表示：
（16）
其中，表示相位频散算子，/>表示振幅补偿算子，/>和/>表示相位校正算子，/>表示高频稳定算子。/>，/>，/>表示模拟频率，表示参考频率，Q表示品质因子，μ表示伴随波场。
求取伴随波场μ可以通过下述公式（17）表示：
（17）
其中，表示伴随震源，/>表示伴随算子。伴随算子/>决定了目标模型的参数的散射方式，不同的波动方程系统对应参数的散射方式也不相同。
针对基于伴随算子和伴随震源，以时间反向的方式求解伴随方程，得到震源对应的逆时伴随波场，可以通过下述公式（18）表示：
（18）
其中，表示相位频散算子，/>表示振幅补偿算子，/>和/>表示相位校正算子，/>表示高频稳定算子。/>，/>，/>表示模拟频率，表示参考频率，Q表示品质因子，μ表示伴随波场，t表示时间。公式（18）为伴随方程。
步骤S104，基于震源波场和逆时伴随波场，确定目标函数的补偿梯度。
在一些实施例中，步骤S104可以包括：将震源波场和逆时伴随波场互相关，得到目标函数的补偿梯度。
针对将震源波场和逆时伴随波场互相关，得到目标函数的补偿梯度G，可以通过下述公式（19）表示：
（19）
其中，C表示目标函数，m表示初始参数，A表示波场延拓算子，μ表示逆时伴随波场，p表示震源波场。
步骤S105，根据补偿梯度对目标函数的初始参数进行迭代更新，直至迭代更新的次数达到预设值。
在一些实施例中，迭代更新的次数对应的预设值可以预先设定，如预设目标函数的初始参数需要迭代n次，n为大于或等于1的正整数。迭代更新可以包括：通过局部优化算法迭代更新初始参数。
在一些实施例中，步骤S105之后，全波形反演方法还可以包括：将最后一次迭代更新后的参数作为目标函数对应的全波形反演结果。
步骤S105可以通过下述公式（20）表示：
（20）
其中，表示第k次迭代更新后的参数，/>表示第k+1次迭代得到的参数，/>表示第k次迭代时的下降步长，/>表示第k次迭代时的下降方向。/>可以由第k次迭代的补偿梯度确定。k为大于或等于0的整数。
作为示例，预先设定迭代更新的次数为3次，首先根据初始参数对应的第一补偿梯度，计算得到第一参数，根据第一参数确定其对应的第二补偿梯度，将初始参数更新为第一参数，将第一补偿梯度更新为第二补偿梯度；根据第二补偿梯度，计算得到第二参数，根据第二参数确定其对应的第三补偿梯度，将第一参数更新为第二参数，将第二补偿梯度更新为第三补偿梯度；根据第三补偿梯度，计算得到第三参数，根据第三参数确定其对应的第四补偿梯度，将第二参数更新为第三参数，将第三补偿梯度更新为第四补偿梯度。并将第三参数作为目标函数对应的全波形反演结果。
在一些实施例中，对所述全波形反演方法的处理流程示意图二，如图2所示，包括：
步骤S201，基于模拟数据，确定对应的第一稳定补偿算子。
步骤S202，基于第一稳定补偿算子，确定模拟数据对应的波场延拓算子。
步骤S203，基于波场延拓算子和震源，以时间正向的方式求解补偿方程，得到震源对应的正传波场。
步骤S204，对正传波场以时间反向的方式进行重构，得到震源波场。
针对步骤S201-S204的每个步骤的具体说明，与上述步骤S102相同，这里不再赘述。
本申请实施例的方法，通过步骤S201-S204，波场延拓算子可解耦表示为不同的独立算子，能够有效解耦表征振幅和相位作用的方程，在振幅补偿的同时不改变相位信息，提高了全波形反演的精度，第一稳定补偿算子中的拉普拉斯算子为固定阶数，因此它自然适应于变化剧烈的衰减介质，减少了全波形反演的计算过程的时间，进而提高了全波形反演的效率。
在一些实施例中，对所述全波形反演方法的处理流程示意图三，如图3所示，包括：
步骤S301，基于模拟数据，确定对应的第二稳定补偿算子。
步骤S302，基于第二稳定补偿算子，求解衰减方程，得到震源对应的粘声波场。
步骤S303，基于粘声波场和观测数据，确定伴随震源。
步骤S304，基于伴随震源，确定逆时伴随波场。
针对步骤S301-S304的每个步骤的具体说明，与上述步骤S103相同，这里不再赘述。
本申请实施例的方法，通过步骤S301-S304，能够有效解耦表征振幅和相位作用的方程，在振幅补偿的同时不改变相位信息，无需引入滤波器，即可稳定补偿振幅，缓解了衰减效应对梯度的影响，衰减和非衰减区的梯度照明得到了有效平衡，进而提高了全波形反演的精度，第二稳定补偿算子中的拉普拉斯算子为固定阶数，因此它自然适应于变化剧烈的衰减介质，减少了全波形反演的计算过程的时间，进而提高了全波形反演的效率。
在一些实施例中，对所述全波形反演方法的处理流程示意图四，如图4所示，包括：
步骤S401，基于伴随震源，确定对应的第三稳定补偿算子。
步骤S402，基于第三稳定补偿算子，确定伴随震源对应的伴随算子。
步骤S403，基于伴随算子和伴随震源，以时间反向的方式求解伴随方程，得到震源对应的逆时伴随波场。
针对步骤S401-S403的每个步骤的具体说明，与上述步骤S103相同，这里不再赘述。
本申请实施例的方法，通过步骤S401-S403，伴随算子可解耦表示为不同的独立算子，能够有效解耦表征振幅和相位作用的方程，在振幅补偿的同时不改变相位信息，提高了全波形反演的精度，第三稳定补偿算子中的拉普拉斯算子为固定阶数，因此它自然适应于变化剧烈的衰减介质，减少了全波形反演的计算过程的时间，进而提高了全波形反演的效率。
在一些实施例中，对所述全波形反演方法的处理流程示意图五，如图5所示，包括：
步骤S501，将震源波场和逆时伴随波场互相关，得到目标函数的补偿梯度。
针对步骤S501的具体说明，与上述步骤S104相同，这里不再赘述。
步骤S502，根据补偿梯度对目标函数的初始参数进行迭代更新，直至迭代更新的次数达到预设值。
步骤S503，将最后一次迭代更新后的参数作为目标函数对应的全波形反演结果。
针对步骤S502-S503的每个步骤的具体说明，与上述步骤S105相同，这里不再赘述。
本申请实施例的方法，通过步骤S501-S503，能够有效解耦表征振幅和相位作用的方程，在振幅补偿的同时不改变相位信息，无需引入滤波器，即可稳定补偿振幅，缓解了衰减效应对梯度的影响，衰减和非衰减区的梯度照明得到了有效平衡，进而提高了全波形反演的精度，减少了全波形反演的计算过程的时间，进而提高了全波形反演的效率。
图6示出了本申请实施例提供的全波形反演方法的一种应用场景图。
参考图6，本申请实施例提供的全波形反演方法的一种应用场景，应用于针对粘声介质的衰减补偿的全波形反演。首先，建立目标函数。设定目标函数的反演参数，反演参数可以为速度参数。沿时间正向解耦补偿方程，获得正传补偿波场，再沿时间反向重构正传补偿波场，得到震源波场；沿时间正向解耦衰减方程，获得粘声波场，并根据粘声波场计算伴随震源，根据伴随震源，沿时间反向解耦伴随方程，得到粘声介质的衰减补偿的逆时伴随波场。将震源波场和逆时伴随波场互相关，计算得到目标函数的补偿梯度。根据补偿梯度对目标函数的反演参数进行迭代更新。判断迭代更新的次数达到最大的迭代次数；若否，则将目标函数的反演参数设定为迭代更新后的反演参数再次进行迭代计算。若是，则将最后一次迭代更新后的参数作为目标函数对应的最终的全波形反演结果。
可以理解，图6的全波形反演方法的应用场景只是本申请实施例中的部分示例性的实施方式，本申请实施例中全波形反演方法的应用场景包括但不限于图6所示的全波形反演方法的应用场景。
图7示出了本申请实施例提供的全波形反演方法的另一种应用场景图。
参考图7，本申请实施例提供的全波形反演方法的另一种应用场景，应用于全波形反演的目标函数下降曲线。
图7的（a）为所有地震记录的全局目标函数下降曲线、（b）为高衰减区域地震记录的目标函数下降曲线、（c）为低减区域地震记录的目标函数下降曲线。其中，虚线表示常规衰减FWI（Full Waveform Inversion, 全波形反演）和补偿FWI。补偿FWI可以包括针对粘声介质的衰减补偿的全波形反演。
通过图7的（a）可知，补偿FWI曲线在前期的迭代中下降速度略快，但是最后与常规FWI下降程度相似。其原因可解释为高衰减区模拟记录与观测记录振幅误差的绝对值低于弱衰减区域，导致低衰减区域数据对误差函数的贡献远大于高衰减区域。因此，误差函数下降曲线将被低衰减区主导，导致实线与虚线最后的目标函数误差曲线接近。
通过对比图7的（b）和图7的（c）可以确定，使用高衰减区域的地震记录进行反演时，补偿FWI的下降速度明显优于常规FWI且误差函数较小，而使用低衰减区域的地震记录时，二者的目标函数下降趋势相近，从而证明了针对粘声介质的衰减补偿的全波形反演可以提高收敛速度，降低时间成本。
可以理解，图7的全波形反演方法的应用场景只是本申请实施例中的部分示例性的实施方式，本申请实施例中全波形反演方法的应用场景包括但不限于图7所示的全波形反演方法的应用场景。
下面继续说明本申请实施例提供的全波形反演装置90的实施为软件模块的示例性结构，在一些实施例中，如图9所示，全波形反演装置90中的软件模块可以包括：设定模块901，可以用于设定震源，并构建震源对应的目标函数；第一确定模块902，可以用于基于震源对应的模拟数据，确定震源对应的震源波场；第二确定模块903，可以用于基于震源对应的模拟数据和观测数据，确定震源对应的逆时伴随波场；梯度确定模块904，可以用于基于震源波场和逆时伴随波场，确定目标函数的补偿梯度；迭代模块905，可以用于根据补偿梯度对目标函数的初始参数进行迭代更新，直至迭代更新的次数达到预设值。
需要说明的是，本申请实施例装置的描述，与上述方法实施例的描述是类似的，具有同方法实施例相似的有益效果，因此不做赘述。对于本申请实施例提供的全波形反演装置中未尽的技术细节，可以根据图1至图7中任一附图的说明而理解。
根据本申请的实施例，本申请还提供了一种电子设备和一种非瞬时计算机可读存储介质。
图9示出了可以用来实施本申请的实施例的示例电子设备800的示意性框图。电子设备旨在表示各种形式的数字计算机，诸如，膝上型计算机、台式计算机、工作台、个人数字助理、服务器、刀片式服务器、大型计算机、和其它适合的计算机。电子设备还可以表示各种形式的移动装置，诸如，个人数字处理、蜂窝电话、智能电话、可穿戴设备和其它类似的计算装置。本文所示的部件、它们的连接和关系、以及它们的功能仅仅作为示例，并且不意在限制本文中描述的和/或者要求的本申请的实现。
如图9所示，电子设备800包括计算单元801，其可以根据存储在只读存储器（ROM）802中的计算机程序或者从存储单元808加载到随机访问存储器（RAM）803中的计算机程序，来执行各种适当的动作和处理。在RAM 803中，还可存储电子设备800操作所需的各种程序和数据。计算单元801、ROM 802以及RAM 803通过总线804彼此相连。输入/输出（I/O）接口805也连接至总线804。
电子设备800中的多个部件连接至I/O接口805，包括：输入单元806，例如键盘、鼠标等；输出单元807，例如各种类型的显示器、扬声器等；存储单元808，例如磁盘、光盘等；以及通信单元809，例如网卡、调制解调器、无线通信收发机等。通信单元809允许电子设备800通过诸如因特网的计算机网络和/或各种电信网络与其他设备交换信息/数据。
计算单元801可以是各种具有处理和计算能力的通用和/或专用处理组件。计算单元801的一些示例包括但不限于中央处理单元（CPU）、图形处理单元（GPU）、各种专用的人工智能（AI）计算芯片、各种运行机器学习模型算法的计算单元、数字信号处理器（DSP）、以及任何适当的处理器、控制器、微控制器等。计算单元801执行上文所描述的各个方法和处理，例如全波形反演方法。例如，在一些实施例中，全波形反演方法可被实现为计算机软件程序，其被有形地包含于机器可读介质，例如存储单元808。在一些实施例中，计算机程序的部分或者全部可以经由ROM 802和/或通信单元809而被载入和/或安装到电子设备800上。当计算机程序加载到RAM 803并由计算单元801执行时，可以执行上文描述的全波形反演方法的一个或多个步骤。备选地，在其他实施例中，计算单元801可以通过其他任何适当的方式（例如，借助于固件）而被配置为执行全波形反演方法。
本文中以上描述的系统和技术的各种实施方式可以在数字电子电路系统、集成电路系统、场可编程门阵列（FPGA）、专用集成电路（ASIC）、专用标准产品（ASSP）、芯片上系统的系统（SOC）、负载可编程逻辑设备（CPLD）、计算机硬件、固件、软件、和/或它们的组合中实现。这些各种实施方式可以包括：实施在一个或者多个计算机程序中，该一个或者多个计算机程序可在包括至少一个可编程处理器的可编程系统上执行和/或解释，该可编程处理器可以是专用或者通用可编程处理器，可以从存储系统、至少一个输入装置、和至少一个输出装置接收数据和指令，并且将数据和指令传输至该存储系统、该至少一个输入装置、和该至少一个输出装置。
用于实施本申请的方法的程序代码可以采用一个或多个编程语言的任何组合来编写。这些程序代码可以提供给通用计算机、专用计算机或其他可编程数据处理装置的处理器或控制器，使得程序代码当由处理器或控制器执行时使流程图和/或框图中所规定的功能/操作被实施。程序代码可以完全在机器上执行、部分地在机器上执行，作为独立软件包部分地在机器上执行且部分地在远程机器上执行或完全在远程机器或服务器上执行。
在本申请的上下文中，机器可读介质可以是有形的介质，其可以包含或存储以供指令执行系统、装置或设备使用或与指令执行系统、装置或设备结合地使用的程序。机器可读介质可以是机器可读信号介质或机器可读储存介质。机器可读介质可以包括但不限于电子的、磁性的、光学的、电磁的、红外的、或半导体系统、装置或设备，或者上述内容的任何合适组合。机器可读存储介质的更具体示例会包括基于一个或多个线的电气连接、便携式计算机盘、硬盘、随机存取存储器（RAM）、只读存储器（ROM）、可擦除可编程只读存储器（EPROM或快闪存储器）、光纤、便捷式紧凑盘只读存储器（CD-ROM）、光学储存设备、磁储存设备、或上述内容的任何合适组合。
为了提供与用户的交互，可以在计算机上实施此处描述的系统和技术，该计算机具有：用于向用户显示信息的显示装置（例如，CRT（阴极射线管）或者LCD（液晶显示器）监视器）；以及键盘和指向装置（例如，鼠标或者轨迹球），用户可以通过该键盘和该指向装置来将输入提供给计算机。其它种类的装置还可以用于提供与用户的交互；例如，提供给用户的反馈可以是任何形式的传感反馈（例如，视觉反馈、听觉反馈、或者触觉反馈）；并且可以用任何形式（包括声输入、语音输入或者、触觉输入）来接收来自用户的输入。
可以将此处描述的系统和技术实施在包括后台部件的计算系统（例如，作为数据服务器）、或者包括中间件部件的计算系统（例如，应用服务器）、或者包括前端部件的计算系统（例如，具有图形用户界面或者网络浏览器的用户计算机，用户可以通过该图形用户界面或者该网络浏览器来与此处描述的系统和技术的实施方式交互）、或者包括这种后台部件、中间件部件、或者前端部件的任何组合的计算系统中。可以通过任何形式或者介质的数字数据通信（例如，通信网络）来将系统的部件相互连接。通信网络的示例包括：局域网（LAN）、广域网（WAN）和互联网。
计算机系统可以包括客户端和服务器。客户端和服务器一般远离彼此并且通常通过通信网络进行交互。通过在相应的计算机上运行并且彼此具有客户端-服务器关系的计算机程序来产生客户端和服务器的关系。服务器可以是云服务器，也可以为分布式系统的服务器，或者是结合了区块链的服务器。
应该理解，可以使用上面所示的各种形式的流程，重新排序、增加或删除步骤。例如，本申请中记载的各步骤可以并行地执行也可以顺序地执行也可以不同的次序执行，只要能够实现本申请公开的技术方案所期望的结果，本文在此不进行限制。
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或隐含地包括至少一个该特征。在本申请的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。
以上所述，仅为本申请的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应以所述权利要求的保护范围为准。
Claims (9)
1.一种全波形反演方法，其特征在于，所述方法包括：
设定震源，并构建所述震源对应的目标函数；
基于所述震源对应的模拟数据，确定所述模拟数据对应的第一稳定补偿算子；所述第一稳定补偿算子包括振幅补偿算子、高频稳定算子、相位频散算子、相位校正算子中的至少之一；
基于所述第一稳定补偿算子，确定所述模拟数据对应的波场延拓算子；
基于所述波场延拓算子和所述震源，以时间正向的方式求解补偿方程，得到所述震源对应的正传波场；
对所述正传波场以时间反向的方式进行重构，得到所述震源对应的震源波场；
基于所述震源对应的模拟数据和观测数据，确定所述震源对应的逆时伴随波场；
基于所述震源波场和所述逆时伴随波场，确定所述目标函数的补偿梯度；
根据所述补偿梯度对所述目标函数的初始参数进行迭代更新，直至所述迭代更新的次数达到预设值。
2.根据权利要求1所述的方法，其特征在于，所述基于所述震源对应的模拟数据和观测数据，确定所述震源对应的逆时伴随波场，包括：
基于所述模拟数据，确定对应的第二稳定补偿算子；所述第二稳定补偿算子包括振幅补偿算子、相位频散算子中的至少之一；
基于所述第二稳定补偿算子，求解衰减方程，得到所述震源对应的粘声波场；
基于所述粘声波场和所述观测数据，确定伴随震源；
基于所述伴随震源，确定所述逆时伴随波场。
3.根据权利要求2所述的方法，其特征在于，所述基于所述伴随震源，确定所述逆时伴随波场，包括：
基于所述伴随震源，确定对应的第三稳定补偿算子；所述第三稳定补偿算子包括振幅补偿算子、高频稳定算子、相位频散算子、相位校正算子中的至少之一；
基于所述第三稳定补偿算子，确定所述伴随震源对应的伴随算子；
基于所述伴随算子和所述伴随震源，以时间反向的方式求解伴随方程，得到所述震源对应的逆时伴随波场。
4.根据权利要求1所述的方法，其特征在于，所述基于所述震源波场和所述逆时伴随波场，确定所述目标函数的补偿梯度，包括：
将所述震源波场和所述逆时伴随波场互相关，得到所述目标函数的补偿梯度。
5.根据权利要求1所述的方法，其特征在于，所述根据所述补偿梯度对所述目标函数的初始参数进行迭代更新，直至所述迭代更新的次数达到预设值之后，所述方法还包括：
将最后一次迭代更新后的参数作为所述目标函数对应的全波形反演结果。
6.根据权利要求1或3所述的方法，其特征在于，所述振幅补偿算子、高频稳定算子、相位频散算子、相位校正算子均包括固定阶数的拉普拉斯算子。
7.一种全波形反演装置，其特征在于，所述全波形反演装置包括：
设定模块，用于设定震源，并构建所述震源对应的目标函数；
第一确定模块，用于基于所述震源对应的模拟数据，确定所述模拟数据对应的第一稳定补偿算子；所述第一稳定补偿算子包括振幅补偿算子、高频稳定算子、相位频散算子、相位校正算子中的至少之一；基于所述第一稳定补偿算子，确定所述模拟数据对应的波场延拓算子；基于所述波场延拓算子和所述震源，以时间正向的方式求解补偿方程，得到所述震源对应的正传波场；对所述正传波场以时间反向的方式进行重构，得到所述震源对应的震源波场；
第二确定模块，用于基于所述震源对应的模拟数据和观测数据，确定所述震源对应的逆时伴随波场；
梯度确定模块，用于基于所述震源波场和所述逆时伴随波场，确定所述目标函数的补偿梯度；
迭代模块，用于根据所述补偿梯度对所述目标函数的初始参数进行迭代更新，直至所述迭代更新的次数达到预设值。
8.一种电子设备，其特征在于，包括：
至少一个处理器；以及
与所述至少一个处理器通信连接的存储器；其中，
所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被所述至少一个处理器执行，以使所述至少一个处理器能够执行权利要求1-6中任一项所述的方法。
9.一种存储有计算机指令的非瞬时计算机可读存储介质，其特征在于，所述计算机指令用于使计算机执行根据权利要求1-6中任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310893372.0A CN116629031A (zh) | 2023-07-19 | 2023-07-19 | 一种全波形反演方法、装置、电子设备及存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310893372.0A CN116629031A (zh) | 2023-07-19 | 2023-07-19 | 一种全波形反演方法、装置、电子设备及存储介质 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116629031A true CN116629031A (zh) | 2023-08-22 |
Family
ID=87638572
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310893372.0A Pending CN116629031A (zh) | 2023-07-19 | 2023-07-19 | 一种全波形反演方法、装置、电子设备及存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116629031A (zh) |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20150362622A1 (en) * | 2014-06-17 | 2015-12-17 | Huseyin Denli | Fast Viscoacoustic and Viscoelastic Full Wavefield Inversion |
| CN105467444A (zh) * | 2015-12-10 | 2016-04-06 | 中国石油天然气集团公司 | 一种弹性波全波形反演方法及装置 |
| CN110579795A (zh) * | 2018-06-08 | 2019-12-17 | 中国海洋大学 | 基于被动源地震波形及其逆时成像的联合速度反演方法 |
| US20200348430A1 (en) * | 2019-05-03 | 2020-11-05 | Exxonmobil Upstream Research Company | Inversion, Migration, and Imaging Related to Isotropic Wave-Mode-Independent Attenuation |
| CN113050160A (zh) * | 2021-03-26 | 2021-06-29 | 中国石油大学(北京) | 一种数据阻尼全波形反演方法、装置和计算机设备 |
| CN113156493A (zh) * | 2021-05-06 | 2021-07-23 | 中国矿业大学 | 一种使用归一化震源的时频域全波形反演方法及装置 |
| CN114966831A (zh) * | 2022-01-28 | 2022-08-30 | 东北石油大学 | 一种基于速度-衰减解耦的粘声全波形反演方法 |
| CN114966837A (zh) * | 2022-05-13 | 2022-08-30 | 中国石油大学(北京) | 基于卷积型w2距离目标函数的全波形反演方法及装置 |
-
2023
- 2023-07-19 CN CN202310893372.0A patent/CN116629031A/zh active Pending
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20150362622A1 (en) * | 2014-06-17 | 2015-12-17 | Huseyin Denli | Fast Viscoacoustic and Viscoelastic Full Wavefield Inversion |
| CN105467444A (zh) * | 2015-12-10 | 2016-04-06 | 中国石油天然气集团公司 | 一种弹性波全波形反演方法及装置 |
| CN110579795A (zh) * | 2018-06-08 | 2019-12-17 | 中国海洋大学 | 基于被动源地震波形及其逆时成像的联合速度反演方法 |
| US20200348430A1 (en) * | 2019-05-03 | 2020-11-05 | Exxonmobil Upstream Research Company | Inversion, Migration, and Imaging Related to Isotropic Wave-Mode-Independent Attenuation |
| CN113050160A (zh) * | 2021-03-26 | 2021-06-29 | 中国石油大学(北京) | 一种数据阻尼全波形反演方法、装置和计算机设备 |
| CN113156493A (zh) * | 2021-05-06 | 2021-07-23 | 中国矿业大学 | 一种使用归一化震源的时频域全波形反演方法及装置 |
| CN114966831A (zh) * | 2022-01-28 | 2022-08-30 | 东北石油大学 | 一种基于速度-衰减解耦的粘声全波形反演方法 |
| CN114966837A (zh) * | 2022-05-13 | 2022-08-30 | 中国石油大学(北京) | 基于卷积型w2距离目标函数的全波形反演方法及装置 |
Non-Patent Citations (1)
| Title |
|---|
| 关珊等: "Q补偿粘滞声波全波形反演", 《第五届油气地球物理学术年会论文集》, pages 1 - 4 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Askan et al. | Full waveform inversion for seismic velocity and anelastic losses in heterogeneous structures | |
| Ma et al. | Modeling of the perfectly matched layer absorbing boundaries and intrinsic attenuation in explicit finite-element methods | |
| US20120257477A1 (en) | Amplitude contrast seismic attribute | |
| Le Bouteiller et al. | An accurate discontinuous Galerkin method for solving point-source Eikonal equation in 2-D heterogeneous anisotropic media | |
| Alkhalifah et al. | An eikonal-based formulation for traveltime perturbation with respect to the source location | |
| CN109946742A (zh) | 一种TTI介质中纯qP波地震数据模拟方法 | |
| Fabien-Ouellet | Seismic modeling and inversion using half-precision floating-point numbers | |
| Wang et al. | Adaptive H∞ Kalman filter based random drift modeling and compensation method for ring laser gyroscope | |
| CN110794458A (zh) | 基于时频分析的含气性检测方法、装置及存储介质 | |
| CN112505775B (zh) | 地震波正演模拟方法、装置、存储介质及处理器 | |
| Gulley et al. | A numerical approach for modelling fault-zone trapped waves | |
| CN116629031A (zh) | 一种全波形反演方法、装置、电子设备及存储介质 | |
| CN109725346B (zh) | 一种时间-空间高阶vti矩形网格有限差分方法及装置 | |
| He et al. | A numerical dispersion-dissipation analysis of discontinuous Galerkin methods based on quadrilateral and triangular elements | |
| WO2022216843A1 (en) | System and method for estimating one-way propagation operators | |
| Mulder | Performance of old and new mass‐lumped triangular finite elements for wavefield modelling | |
| Roux et al. | Stability/dispersion analysis of the discontinuous Galerkin linearized shallow‐water system | |
| CN114002741B (zh) | 叠前深度偏移方法及装置、计算机可读存储介质 | |
| CN117849875B (zh) | 一种地震信号分析方法、系统、装置及存储介质 | |
| CN112241022B (zh) | 基于射线密度生成层析反演模型速度界面的方法和装置 | |
| CN112114364B (zh) | 偶极横波反射波补偿方法及装置 | |
| Jing et al. | An optimized time-space-domain finite difference method with piecewise constant interpolation coefficients for scalar wave propagation | |
| CN112859170B (zh) | 一种基于时域高阶有限差分法的高精度波场数值模拟方法 | |
| CN112782756B (zh) | 一种基于自适应构造约束的约束层速度反演方法及系统 | |
| CN111948706B (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 |