CN110858002B - 拓展显式差分稳定性条件的波场模拟方法、装置及设备 - Google Patents

拓展显式差分稳定性条件的波场模拟方法、装置及设备 Download PDF

Info

Publication number
CN110858002B
CN110858002B CN201810967522.7A CN201810967522A CN110858002B CN 110858002 B CN110858002 B CN 110858002B CN 201810967522 A CN201810967522 A CN 201810967522A CN 110858002 B CN110858002 B CN 110858002B
Authority
CN
China
Prior art keywords
wave field
new
eigenvalue
matrix
wave
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
CN201810967522.7A
Other languages
English (en)
Other versions
CN110858002A (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.)
Institute of Geology and Geophysics of CAS
Original Assignee
Institute of Geology and Geophysics of CAS
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 Institute of Geology and Geophysics of CAS filed Critical Institute of Geology and Geophysics of CAS
Priority to CN201810967522.7A priority Critical patent/CN110858002B/zh
Priority to PCT/CN2019/083710 priority patent/WO2020038007A1/zh
Publication of CN110858002A publication Critical patent/CN110858002A/zh
Application granted granted Critical
Publication of CN110858002B publication Critical patent/CN110858002B/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/30—Analysis
    • G01V1/308—Time lapse or 4D effects, e.g. production related effects to the formation
    • 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/32—Transforming one recording into another or one representation into another

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Fluid Mechanics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及一种拓展显式差分稳定性条件的波场模拟方法、装置及设备。所述方法包括：根据设定的时间步长和波场模拟模型参数对波动方程进行有限差分离散，得到离散后的波动方程；对所述离散后的波动方程的震源函数进行时间频散变换得到新的震源函数；对所述离散后的波动方程的增长矩阵进行特征值扰动得到新的增长矩阵；将所述新的震源函数和所述新的增长矩阵代入所述离散后的波动方程，得到新的稳定的波场迭代式；采用所述新的波场迭代式进行波场数值模拟得到波场数据；针对所述波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据；存储所述的逆时间频散变换处理之后的波场数据。

Description

拓展显式差分稳定性条件的波场模拟方法、装置及设备
技术领域
本发明涉及波场模拟领域，具体涉及拓展显式差分稳定性条件的波场模拟方法、装置及设备。
背景技术
显式有限差分由于其形式简单易于实现，是目前波场模拟领域应用最广泛的时间偏导数求解方法。其中模拟精度和计算效率是显式有限差分所关心的最重要的两个问题。而时间步长的选取会同时影响数值模拟的精度和效率。
足够小的时间步长即使在长时程数值模拟中仍然可以得到高精度的数值模拟结果而不出现时间频散，但是小的时间步长意味着更多的计算量，计算效率较低；与之相反，显式差分的最大时间步长受到Courant-Friedrichs-Lewy(CFL)稳定性条件的限制，采用接近于CFL稳定性条件上限的大时间步长时，计算效率较高，但是计算精度较低，即使在短时程的模拟结果中仍然会出现较为严重时间频散。于是，在具有精细结构、高速体的速度模型中，由于只能采用小时间步长而增加了数值模拟的迭代次数，降低了计算效率。
时间频散被证明独立于空间频散，人们提出了多种去除时间频散的方法。这使得我们可以采用更大的时间步长(接近于稳定条件上限)进行数值模拟而不用担心时间频散。虽然大时间步长的时间频散问题得到了解决，截至目前，波场模拟领域显式有限差分的时间步长仍然受到稳定条件的限制而不能取得过大。
隐式有限差分是无条件稳定的，不考虑计算精度的情况下，可以采用任意的时间步长而不受CFL稳定性条件的限制。隐式有限差分法需要求解大型的矩阵，比显示差分需要更多的计算量，这类方法的计算效率较低。
电磁波模拟领域出现了采用特征值扰动的思路实现显式有限差分无条件稳定的方法，这种数值模拟方法还没有应用于波场模拟领域。
发明内容
本发明提出拓展显式差分稳定性条件的波场模拟方法、装置及设备，旨在波场模拟领域采用特征值扰动的方法，实现拓展显式有限差分稳定性条件，并且将时间频散变换和逆时间频散变换引入以去除大时间步长模拟所产生的时间频散，以提供一种高精度、高效率的波场数值模拟新方法、新工具。
本发明提供一种拓展显式差分稳定性条件的波场模拟方法，包括：
根据设定的时间步长和波场模拟模型参数对波动方程进行有限差分离散，得到离散后的波动方程；
对所述离散后的波动方程的震源函数进行时间频散变换得到新的震源函数；
对所述离散后的波动方程的增长矩阵进行特征值扰动得到新的增长矩阵；
将所述新的震源函数和所述新的增长矩阵代入所述离散后的波动方程，得到新的稳定的波场迭代式；
采用所述新的波场迭代式进行波场数值模拟得到波场数据；
针对所述波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据；
存储所述的逆时间频散变换处理之后的波场数据。
进一步地，所述波场模拟模型参数包括：空间网格大小，模型全局速度、时间步长、震源函数、时间和空间有限差分的离散格式。
进一步地，所述对所述离散后的波动方程的震源函数进行时间频散变换得到新的震源函数，包括：
对所述离散后的波动方程的震源函数利用时间频散变换关系计算理论相移；
基于所述理论相移对所述震源函数进行变换，该变换可以看做修正的傅里叶变换；
继续进行逆傅里叶变换，得到新的震源函数。
进一步地，所述对所述离散后的波动方程的增长矩阵进行特征值扰动得到新的增长矩阵，包括：
对所述波动方程的增长矩阵进行特征值分解，得到所述增长矩阵的特征向量矩阵和特征值对角矩阵；所述特征值对角矩阵由所述增长矩阵的特征值组成；
采用特征值扰动的方法对所述增长矩阵中模大于1的特征值进行扰动，使这部分特征值的模变为1，得到新特征值；
将所述增长矩阵中模小于等于1的特征值和所述新特征值组成新的对角矩阵；
采用所述新的对角矩阵和原来的所述特征向量矩阵构建新的增长矩阵。
进一步地，所述针对所述波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据，包括：
利用时间离散之后的频散变换关系计算所述目的波场数据的实际相移；
基于所述实际相移对所述目的波场数据进行变换，该变换可以看做修正的傅里叶变换；
继续进行逆傅里叶变换，得到消除时间频散的波场数据。
本发明实施例还提供一种拓展显式差分稳定性条件的波场模拟装置，其特征在于，所述装置包括：
离散单元，根据设定的时间步长和波场模拟模型参数对波动方程进行有限差分离散，得到离散后的波动方程；
震源函数变换单元，对所述离散后的波动方程的震源函数进行时间频散变换得到新的震源函数；
特征值扰动单元，对所述离散后的波动方程的增长矩阵进行特征值扰动得到新的增长矩阵；
迭代式构成单元，所述新的震源函数和所述新的增长矩阵代入所述离散后的波动方程，得到新的稳定的波场迭代式；
波场模拟单元，采用所述新的波场迭代式进行波场数值模拟得到波场数据；
逆时间频散变换单元，针对所述波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据；
存储单元，存储所述的逆时间频散变换处理之后的波场数据。
进一步，所述震源函数变换单元包括：
第一计算模块，对所述离散后的波动方程的震源函数利用时间频散变换关系计算理论相移；
第一修正的傅里叶变换模块，基于所述理论相移对所述震源函数进行变换，该变换可以看做修正的傅里叶变换；
第一逆傅里叶变换模块，继续进行逆傅里叶变换，得到新的震源函数；
所述逆时间频散变换单元包括：
第二计算模块，利用时间频散变换关系计算目的波场数据的实际相移；
第二修正的傅里叶变换模块，基于所述实际相移对所述目的波场数据进行变换，该变换可以看做修正的傅里叶变换；
第二逆傅里叶变换模块，继续进行逆傅里叶变换，得到新的目的波场数据。
进一步地，所述特征值扰动单元包括：
特征值分解模块，对所述波动方程的增长矩阵进行特征值分解，得到所述增长矩阵的特征向量矩阵和特征值对角矩阵；所述特征值对角矩阵由所述增长矩阵的特征值组成；
扰动模块，采用特征值扰动的方法对所述增长矩阵中模大于1的特征值进行扰动，使这部分特征值的模变为1，得到新特征值；
对角矩阵构成模块，将所述增长矩阵中模小于等于1的特征值和所述新特征值组成新的对角矩阵；
增长矩阵构成模块，采用所述新的对角矩阵和原来的所述特征向量矩阵构建新的增长矩阵。
本发明还提供一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现上述所述的方法。
本发明还提供一种计算机可读存储介质，其上存储有处理器程序,其特征在于，该处理器程序用于执行上述所述拓展显式差分稳定性条件的波场模拟方法。
本发明实施例采用特征值扰动的思路，针对地震波场数值模拟发展了一种显式无条件稳定的模拟方法，极大地提高了数值模拟的效率；同时引入时间频散变换和逆时间频散变换来去除时间频散，同时保证了数值模拟效率和模拟精度。
附图说明
图1是本发明一实施例提供的一种拓展显式差分稳定性条件的波场模拟方法流程示意图；
图2是本发明一实施例提供的一种震源函数时间频散变换前后对比图；
图3是本发明一实施例提供的特征值扰动示意图；
图4是本发明一实施例提供的不同时间步长的单道波形记录对比示意图；
图5是本发明一实施例提供的一种拓展显式差分稳定性条件的波场模拟装置组成示意图；
图6是本发明另一实施例提供的一种拓展显式差分稳定性条件的波场模拟装置组成示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，以下将结合附图和实施例，对本发明技术方案的具体实施方式进行更加详细、清楚的说明。然而，以下描述的具体实施方式和实施例仅是说明的目的，而不是对本发明的限制。其只是包含了本发明一部分实施例，而不是全部的实施例，本领域技术人员对于本发明的各种变化获得的其他实施例，都属于本发明保护的范围。
应该理解的是，虽然第一、第二、第三等用语可使用于本文中用来描述各种元件或组件，但这些元件或组件不应被这些用语所限制。这些用语仅用以区分一个元件或组件与另一元件或组件。因此，下述讨论之第一元件或组件，在不脱离本发明之内容下，可被称为第二元件或第二组件。
CFL条件是以Courant，Friedrichs，Lewy三个人的名字命名的，他们最早在1928年一篇关于偏微分方程的有限差分方法的文章中首次提出这个概念的时候，并不是用来分析差分格式的稳定性，而是仅仅以有限差分方法作为分析工具来证明某些偏微分方程的解的存在性的。其基本思想是先构造差分方程得到一个逼近解的序列，只要知道在给定的网格系统下这个逼近序列收敛，那么就很容易证明这个收敛解就是原微分方程的解。
要使这个逼近序列收敛，必须满足一个条件，那就是著名的CFL条件，记述如下：一个数值方法只有在其依赖的数值域内包含双曲型方程的差分格式收敛的必要条件是差分格式的依赖域包含了微分方程的依赖域。
随着计算机的迅猛发展，有限差分方法和有限体积方法越来越多的应用于流体力学的数值模拟中，CFL条件作为一个格式稳定性和收敛性的判据，也随之显得非常重要了。但值得注意的是，CFL条件仅仅是稳定性(收敛性)的必要条件，而不是充分条件。
图1是本发明一实施例提供的一种拓展显式差分稳定性条件的波场模拟方法流程示意图，包括以下步骤。
在步骤S110中，根据设定的时间步长和波场模拟模型参数对波动方程进行有限差分离散，得到离散后的波动方程。
波场模拟模型参数包括：空间网格大小，模型全局速度、时间步长、震源函数、时间和空间有限差分的离散格式。
以二维标量波动方程为例：
Figure BDA0001775337080000071
其中u(x,z,t)是波场值，c(x,z)是波场传播速度，s(xs,zs,t)是震源函数。
时间偏导采用2阶有限差分进行离散，空间偏导采用4阶有限差分进行离散。离散后的波动方程可以写为矩阵形式：
Figure BDA0001775337080000072
其中Un+1和Un分别代表(n+1)Δt和nΔt时刻的波场值列向量；Vn+1和Vn是引入的辅助变量，分别表示nΔt和(n-1)Δt时刻的波场值列向量，即Vn+1＝Un，Vn＝Un-1。I是单位对角矩阵，A是增长矩阵。Sn是包含震源时间函数的空间位置的列向量。
在步骤S120中，对离散后的波动方程的震源函数进行时间频散变换得到新的震源函数。
图2是本发明一实施例提供的一种震源函数时间频散变换前后对比图。震源采用主频为20Hz的雷克子波，时间步长采用7ms。图中曲线121是原始震源函数，经过时间频散变换操作之后，可以得到图中的曲线122，该曲线中被加入了时间频散，以便于模拟之后得到的波场中的时间频散可以通过逆时间频散变换得以准确地滤除。
采用大时间步长进行波场模拟会产生时间频散，降低数值模拟的精度。时间频散变换和逆时间频散变换分别一种预处理和后处理的技术，波场模拟开始之前对震源函数进行时间频散变换加入预测的时间频散，波场模拟完成之后对地震道采用逆时间频散变换去除时间频散，可以显著地提高模拟精度。以时间2阶有限差分为例，该步骤包括以下子步骤S121、S122、S123。
步骤S121，对离散后的波动方程的震源函数利用时间离散之后的频散关系计算理论相移。
计算公式为
Figure BDA0001775337080000081
其中，ω是有效频率，Δt是时间步长。
步骤S122，基于理论相移对震源函数进行变换，该变换可以看做修正的傅里叶变换。
具体的变换公式为：
Figure BDA0001775337080000082
其中S(t)是震源函数。
步骤S123，继续进行逆傅里叶变换，得到新的震源函数。
具体的变换公式为：
Figure BDA0001775337080000083
在步骤S130中，对离散后的波动方程的增长矩阵进行特征值扰动得到新的增长矩阵。
数学分析和实验表明，给定一个在CFL稳定性条件限定内的时间步长，对应增长矩阵的特征值都分布在单位圆上；如果时间步长超过了CFL稳定性条件的限制，会有部分特征值分布在单位圆之外，这部分特征值会导致不稳定现象。特征值扰动可以对不稳定的特征值进行扰动，从而保证采用任意时间步长时的增长矩阵的稳定性。
针对增长矩阵特征值的扰动可以在波场模拟开始之前进行，采用特征值扰动之后的增长矩阵进行波场模拟就不会产生不稳定现象。特征值扰动法作为一个预处理的过程，能够大幅提高波场模拟的效率。与此同时，有效波场分量的频带宽度主要分布在不需要被扰动的稳定特征值所对应的频带范围，因此特征值扰动并不会影响有效波场。该步骤包括以下子步骤S131、S132、S133、S134。
步骤S131，对波动方程的增长矩阵进行特征值分解，得到增长矩阵的特征向量矩阵U和特征值对角矩阵；特征值对角矩阵由增长矩阵的特征值组成。
对矩阵A进行特征值分解：A＝UΛU-1，其中特征向量矩阵U包含了增长矩阵A的特征向量，Λ是由矩阵A的特征值组成的特征值对角矩阵。
步骤S132，采用特征值扰动的方法对增长矩阵中模大于1的特征值进行扰动，使这部分特征值的模变为1，得到新特征值。
扰动公式为：
Figure BDA0001775337080000091
其中扰动系数γ＝1。λi为增长矩阵A中模大于1的特征值。
图3是本发明一实施例提供的特征值扰动示意图。如图3所示，a图、b图中采用时间步长0.5ms和1.5ms时，均满足稳定性条件，可以看到特征值以圆圈表示，分布在单位圆上。c图中，时间步长采用2ms时，由于给定的模型参数稳定性上限是Δt＝1.530ms，那么2ms的时间步长超过了稳定性条件，可以看到有部分特征值分布在单位圆外面。通过特征值扰动，将不稳定特征值的模扰动为1可以得到d图。
步骤S133，将增长矩阵中模小于等于1的特征值和新特征值组成新的对角矩阵
Figure BDA0001775337080000101
步骤S134，采用新的对角矩阵
Figure BDA0001775337080000102
和原来的特征向量矩阵U构建新的增长矩阵
Figure BDA0001775337080000103
构建新的增长矩阵
Figure BDA0001775337080000104
在步骤S140中，将新的震源函数和新的增长矩阵代入离散后的波动方程，得到新的稳定的波场迭代式。
特征值扰动之后得到新的稳定的波场迭代式：
Figure BDA0001775337080000105
其中
Figure BDA0001775337080000106
中的震源时间函数是时间频散变换之后的S′(t)。
在步骤S150中，采用新的波场迭代式进行波场数值模拟得到波场数据。
在步骤S160中，针对波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据。该步骤包括字步骤S161、S162、S163。
步骤S161，利用时间离散之后的频散变换关系计算目的波场数据的实际相移。
计算公式为：
Figure BDA0001775337080000107
其中，θ(ω,Δt)是实际相移，ω是有效频率，Δt是时间步长。
步骤S162，基于实际相移对目的波场数据进行变换，该变换可以看做修正的傅里叶变换。
变换公式为：
Figure BDA0001775337080000108
其中u(t)是需要输出的目的波场数据。
步骤S163，继续进行逆傅里叶变换，得到消除时间频散的波场数据u′(t)。
变换公式为：
Figure BDA0001775337080000111
在步骤S170中，存储的逆时间频散变换处理之后的波场数据。
图4是本发明一实施例提供的不同时间步长的单道波形记录对比示意图。如图4所示，a图中虚线是理论波场，实线是经过震源函数时间频散变换和特征值扰动处理的波场模拟后的波形。分别采用了时间步长Δt＝2，3，4，5，6和7ms进行波场模拟。b图是对a图的波形进行逆时间频散变换处理之后的波形。从b图可知，经过上述方法处理之后的实际波形与理论波形基本重合，取得了良好的恢复效果。7ms时已超过奈奎斯特(Nyquist)采样上限，所以没有得到恢复的波形。
由于显式有限差分的最大时间步长受到CFL稳定性条件的限制，在本实施例中，给定的模型参数稳定性上限是Δt＝1.530ms，而本发明实施例提供的方法可以突破稳定性条件直到Nyquist采样上限。
图5是本发明一实施例提供的一种拓展显式差分稳定性条件的波场模拟装置组成示意图。所述装置包括离散单元1、震源函数变换单元2、特征值扰动单元3、迭代式构成单元4、波场模拟单元5、逆时间频散变换单元6、存储单元7。
离散单元1根据设定的时间步长和波场模拟模型参数对波动方程进行有限差分离散，得到离散后的波动方程。震源函数变换单元2对离散后的波动方程的震源函数进行时间频散变换得到新的震源函数。特征值扰动单元3对离散后的波动方程的增长矩阵进行特征值扰动得到新的增长矩阵。迭代式构成单元4将新的震源函数和新的增长矩阵代入离散后的波动方程，得到新的稳定的波场迭代式。波场模拟单元5采用新的波场迭代式进行波场数值模拟得到波场数据。逆时间频散变换单元6针对波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据。存储单元7存储的逆时间频散变换处理之后的波场数据。
图6是本发明另一实施例提供的一种拓展显式差分稳定性条件的波场模拟装置组成示意图。
在本实施例中，震源函数变换单元2包括第一计算模块21、第一修正的傅里叶变换模块22、第一逆傅里叶变换模块23。
第一计算模块21对离散后的波动方程的震源函数利用时间频散变换关系计算理论相移。第一修正的傅里叶变换模块22基于理论相移对震源函数进行修正的傅里叶变换。第一逆傅里叶变换模块23继续对傅里叶变换后的震源函数进行逆傅里叶变换，得到新的震源函数。
特征值扰动单元3包括特征值分解模块31、扰动模块32、对角矩阵构成模块33、增长矩阵构成模块34。
特征值分解模块31对波动方程的增长矩阵进行特征值分解，得到增长矩阵的特征向量矩阵和特征值对角矩阵；特征值对角矩阵由增长矩阵的特征值组成。扰动模块32采用特征值扰动的方法对增长矩阵中模大于1的特征值进行扰动，使这部分特征值的模变为1，得到新特征值。对角矩阵构成模块33将增长矩阵中模小于等于1的特征值和新特征值组成新的对角矩阵。增长矩阵构成模块34采用新的对角矩阵和原来的特征向量矩阵构建新的增长矩阵。
逆时间频散变换单元6包括第二计算模块61、第二修正的傅里叶变换模块62、第二逆傅里叶变换模块63。
第二计算模块61利用时间频散变换关系计算目的波场数据的实际相移。第二修正的傅里叶变换模块62基于实际相移对目的波场数据进行修正的傅里叶变换。第二逆傅里叶变换模块63对傅里叶变换后的目的波场数据继续进行逆傅里叶变换，得到新的目的波场数据。
一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，处理器执行程序时实现如上的拓展显式差分稳定性条件的波场模拟方法。
一种计算机可读存储介质，其上存储有处理器程序,该处理器程序用于执行上述的拓展显式差分稳定性条件的波场模拟方法。
其中，计算机可读存储介质可以是只读存储器(ROM)、随机存取存储器(RAM)、光盘只读存储器(CD-ROM)、磁带、软盘和光数据存储设备等，根据实际情况选择，并不以此为限。
需要说明的是，本发明提供的方法和装置，虽然是基于时间显式有限差分和空间有限差分说明的，但是实际应用中，还可以用于伪谱法。
需要说明的是，以上参照附图所描述的各个实施例仅用以说明本发明而非限制本发明的范围，本领域的普通技术人员应当理解，在不脱离本发明的精神和范围的前提下对本发明进行的修改或者等同替换，均应涵盖在本发明的范围之内。此外，除上下文另有所指外，以单数形式出现的词包括复数形式，反之亦然。另外，除非特别说明，那么任何实施例的全部或一部分可结合任何其它实施例的全部或一部分来使用。

Claims (10)

1.一种拓展显式差分稳定性条件的波场模拟方法，包括：
根据设定的时间步长和波场模拟模型参数对波动方程进行有限差分离散，得到离散后的波动方程；
对所述离散后的波动方程的震源函数进行时间频散变换得到新的震源函数；
对所述波动方程的增长矩阵进行特征值分解，获取所述增长矩阵的特征向量矩阵和由所述增长矩阵的特征值组成的特征值对角矩阵并对所述离散后的波动方程的增长矩阵进行特征值扰动，以生成新的对角矩阵，并基于所述新的对角矩阵和原来的所述特征向量矩阵得到新的增长矩阵；
将所述新的震源函数和所述新的增长矩阵代入所述离散后的波动方程，得到新的稳定的波场迭代式；
采用所述新的波场迭代式进行波场数值模拟得到波场数据；
针对所述波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据；
存储所述的逆时间频散变换处理之后的波场数据。
2.根据权利要求1所述的方法，其特征在于，所述波场模拟模型参数包括：空间网格大小，模型全局速度、时间步长、震源函数、时间和空间有限差分的离散格式。
3.根据权利要求1所述的方法，其特征在于，所述对所述离散后的波动方程的震源函数进行时间频散变换得到新的震源函数，包括：
对所述离散后的波动方程的震源函数利用时间频散变换关系计算理论相移；
基于所述理论相移对所述震源函数进行变换，该变换可以看做修正的傅里叶变换；
继续进行逆傅里叶变换，得到新的震源函数。
4.根据权利要求1所述的方法，其特征在于，所述对所述离散后的波动方程的增长矩阵进行特征值扰动，包括：
采用特征值扰动的方法对所述增长矩阵中模大于1的特征值进行扰动，使这部分特征值的模变为1，得到新特征值；
将所述增长矩阵中模小于等于1的特征值和所述新特征值组成所述新的对角矩阵。
5.根据权利要求1所述的方法，其特征在于，所述针对所述波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据，包括：
利用时间离散之后的频散变换关系计算所述目的波场数据的实际相移；
基于所述实际相移对所述目的波场数据进行变换，该变换可以看做修正的傅里叶变换；
继续进行逆傅里叶变换，得到消除时间频散的波场数据。
6.一种拓展显式差分稳定性条件的波场模拟装置，其特征在于，所述装置包括：
离散单元，根据设定的时间步长和波场模拟模型参数对波动方程进行有限差分离散，得到离散后的波动方程；
震源函数变换单元，对所述离散后的波动方程的震源函数进行时间频散变换得到新的震源函数；
特征值扰动单元，对所述波动方程的增长矩阵进行特征值分解，获取所述增长矩阵的特征向量矩阵和由所述增长矩阵的特征值组成的特征值对角矩阵并对所述离散后的波动方程的增长矩阵进行特征值扰动，以生成新的对角矩阵，并基于所述新的对角矩阵和原来的所述特征向量矩阵得到新的增长矩阵；
迭代式构成单元，所述新的震源函数和所述新的增长矩阵代入所述离散后的波动方程，得到新的稳定的波场迭代式；
波场模拟单元，采用所述新的波场迭代式进行波场数值模拟得到波场数据；
逆时间频散变换单元，针对所述波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据；
存储单元，存储所述的逆时间频散变换处理之后的波场数据。
7.根据权利要求6所述的装置，其特征在于，所述震源函数变换单元包括：
第一计算模块，对所述离散后的波动方程的震源函数利用时间频散变换关系计算理论相移；
第一修正的傅里叶变换模块，基于所述理论相移对所述震源函数进行变换，该变换可以看做修正的傅里叶变换；
第一逆傅里叶变换模块，继续进行逆傅里叶变换，得到新的震源函数；
所述逆时间频散变换单元包括：
第二计算模块，利用时间频散变换关系计算目的波场数据的实际相移；
第二修正的傅里叶变换模块，基于所述实际相移对所述目的波场数据进行变换，该变换可以看做修正的傅里叶变换；
第二逆傅里叶变换模块，继续进行逆傅里叶变换，得到新的目的波场数据。
8.根据权利要求6所述的装置，其特征在于，所述特征值扰动单元包括：
特征值分解模块，对所述波动方程的增长矩阵进行特征值分解，得到所述增长矩阵的特征向量矩阵和所述特征值对角矩阵；所述特征值对角矩阵由所述增长矩阵的特征值组成；
扰动模块，采用特征值扰动的方法对所述增长矩阵中模大于1的特征值进行扰动，使这部分特征值的模变为1，得到新特征值；
对角矩阵构成模块，将所述增长矩阵中模小于等于1的特征值和所述新特征值组成所述新的对角矩阵；
增长矩阵构成模块，采用所述新的对角矩阵和原来的所述特征向量矩阵构建新的增长矩阵。
9.一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现如权利要求1至5之任一项所述的方法。
10.一种计算机可读存储介质，其上存储有处理器程序,其特征在于，该处理器程序用于执行上述权利要求1-5任一项所述拓展显式差分稳定性条件的波场模拟方法。
CN201810967522.7A 2018-08-23 2018-08-23 拓展显式差分稳定性条件的波场模拟方法、装置及设备 Active CN110858002B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN201810967522.7A CN110858002B (zh) 2018-08-23 2018-08-23 拓展显式差分稳定性条件的波场模拟方法、装置及设备
PCT/CN2019/083710 WO2020038007A1 (zh) 2018-08-23 2019-04-22 拓展显式差分稳定性条件的波场模拟方法、装置及设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810967522.7A CN110858002B (zh) 2018-08-23 2018-08-23 拓展显式差分稳定性条件的波场模拟方法、装置及设备

Publications (2)

Publication Number Publication Date
CN110858002A CN110858002A (zh) 2020-03-03
CN110858002B true CN110858002B (zh) 2022-07-15

Family

ID=69592266

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810967522.7A Active CN110858002B (zh) 2018-08-23 2018-08-23 拓展显式差分稳定性条件的波场模拟方法、装置及设备

Country Status (2)

Country Link
CN (1) CN110858002B (zh)
WO (1) WO2020038007A1 (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112214938B (zh) * 2020-09-16 2024-04-09 北京动力机械研究所 一种基于dmd方法的流动分离控制系统及方法
CN112505775B (zh) * 2020-10-16 2022-07-26 中国石油大学(华东) 地震波正演模拟方法、装置、存储介质及处理器
CN113742644B (zh) * 2021-09-15 2023-07-04 南方科技大学 基于带宽矩阵的波场值求解的方法、装置及相关设备
CN114357831B (zh) * 2021-12-29 2023-01-24 中海石油(中国)有限公司 无网格广义有限差分正演方法、装置、存储介质及设备
CN116738128A (zh) * 2022-03-04 2023-09-12 本源量子计算科技（合肥）股份有限公司 一种利用量子线路求解含时偏微分方程的方法及装置
CN115373020B (zh) * 2022-08-22 2023-09-29 吉林大学 一种基于离散小波矩量法的地震散射波场数值模拟方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103699798A (zh) * 2013-12-25 2014-04-02 中国石油天然气股份有限公司 一种实现地震波场数值模拟方法
CN105911584A (zh) * 2015-09-25 2016-08-31 中国科学院地质与地球物理研究所 一种隐式交错网格有限差分弹性波数值模拟方法及装置
CN107942375A (zh) * 2017-11-17 2018-04-20 河海大学 一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法
CN107976710A (zh) * 2017-11-17 2018-05-01 河海大学 一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法
GB201805986D0 (en) * 2018-04-11 2018-05-23 Statoil Asa Methods and systems for finite-difference wave equation modelling

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9910174B2 (en) * 2014-07-25 2018-03-06 Seoul National University R&Db Foundation Seismic imaging apparatus and method for performing iterative application of direct waveform inversion
CN106201644A (zh) * 2015-04-29 2016-12-07 中国石油化工股份有限公司 基于gpu的地震数值模拟方法、装置和系统
CN106291664B (zh) * 2015-05-22 2018-11-20 中国石油化工股份有限公司 盲源地震波场模拟方法及系统
CN105044771B (zh) * 2015-08-05 2017-10-27 北京多分量地震技术研究院 基于有限差分法的三维tti双相介质地震波场数值模拟方法
US10621266B2 (en) * 2015-08-25 2020-04-14 Saudi Arabian Oil Company Three-dimensional elastic frequency-domain iterative solver for full waveform inversion
CN105425298B (zh) * 2015-11-11 2018-05-04 中国石油天然气集团公司 一种消除有限差分正演过程中数值频散的方法和装置
CN106709191A (zh) * 2016-12-29 2017-05-24 中国石油大学(北京) 一种地震波场数值模拟方法及装置

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103699798A (zh) * 2013-12-25 2014-04-02 中国石油天然气股份有限公司 一种实现地震波场数值模拟方法
CN105911584A (zh) * 2015-09-25 2016-08-31 中国科学院地质与地球物理研究所 一种隐式交错网格有限差分弹性波数值模拟方法及装置
CN107942375A (zh) * 2017-11-17 2018-04-20 河海大学 一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法
CN107976710A (zh) * 2017-11-17 2018-05-01 河海大学 一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法
GB201805986D0 (en) * 2018-04-11 2018-05-23 Statoil Asa Methods and systems for finite-difference wave equation modelling

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于逆时间频散变换和三阶辛积分的长时程数值模拟;高英杰等;《中国地球科学联合学术年会2016》;20161231;第1190-1192页 *

Also Published As

Publication number Publication date
WO2020038007A1 (zh) 2020-02-27
CN110858002A (zh) 2020-03-03

Similar Documents

Publication Publication Date Title
CN110858002B (zh) 拓展显式差分稳定性条件的波场模拟方法、装置及设备
Tóth et al. On the state-space realization of LPV input-output models: Practical approaches
US10054714B2 (en) Fast viscoacoustic and viscoelastic full wavefield inversion
Deniz et al. Revisiting four approximation methods for fractional order transfer function implementations: Stability preservation, time and frequency response matching analyses
Pettersson et al. A stochastic Galerkin method for the Euler equations with Roe variable transformation
Shivanian Analysis of meshless local and spectral meshless radial point interpolation (MLRPI and SMRPI) on 3-D nonlinear wave equations
Olsen-Kettle Numerical solution of partial differential equations
Roberts Two-dimensional analysis of an iterative nonlinear optimal control algorithm
Liu et al. The almost sure asymptotic stability and $ p $ th moment asymptotic stability of nonlinear stochastic differential systems with polynomial growth
CN112805598A (zh) 拓展有限差分稳定性条件的波场模拟方法、设备及介质
CN108828659B (zh) 基于傅里叶有限差分低秩分解的地震波场延拓方法及装置
Imran et al. Frequency weighted model order reduction technique and error bounds for discrete time systems
Abreu et al. Spectral discretizations analysis with time strong stability preserving properties for pseudo-parabolic models
Cai et al. Finite Frequency Analysis and Synthesis for Singularly Perturbed Systems
Lu et al. Finite element method of BBM-Burgers equation with dissipative term based on adaptive moving mesh
Baffet et al. Adaptive spectral decompositions for inverse medium problems
Khoromskaia et al. Prospects of tensor-based numerical modeling of the collective electrostatics in many-particle systems
Ghaffar et al. Multigrid method for solution of 3D Helmholtz equation based on HOC schemes
Kanaun et al. Scattering of plane monochromatic waves from a heterogeneous inclusion of arbitrary shape in a poroelastic medium: An efficient numerical solution
Saibaba et al. Fast Kalman filter using hierarchical matrices and a low-rank perturbative approach
Taunay et al. Quadrature-based moment methods for kinetic plasma simulations
Dirckx et al. Frequency extraction for bem matrices arising from the 3d scalar helmholtz equation
Axelsson Milestones in the development of iterative solution methods
CN112630823B (zh) 基于交错网格低秩有限差分的三维弹性波场数值模拟方法及系统
Tian et al. Exponential stability of switched positive homogeneous systems

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