CN106842306B - 一种全局优化的交错网格有限差分正演模拟方法和装置 - Google Patents
一种全局优化的交错网格有限差分正演模拟方法和装置 Download PDFInfo
- Publication number
- CN106842306B CN106842306B CN201710253444.XA CN201710253444A CN106842306B CN 106842306 B CN106842306 B CN 106842306B CN 201710253444 A CN201710253444 A CN 201710253444A CN 106842306 B CN106842306 B CN 106842306B
- Authority
- CN
- China
- Prior art keywords
- difference
- finite
- mesh finite
- staggering mesh
- staggering
- 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.)
- Expired - Fee Related
Links
- 238000005457 optimization Methods 0.000 title claims abstract description 125
- 238000000034 method Methods 0.000 title claims abstract description 64
- 238000001228 spectrum Methods 0.000 claims abstract description 77
- 238000004088 simulation Methods 0.000 claims abstract description 56
- 230000006870 function Effects 0.000 claims description 144
- 238000002922 simulated annealing Methods 0.000 claims description 54
- 238000005070 sampling Methods 0.000 claims description 25
- 238000009795 derivation Methods 0.000 claims description 12
- 238000004422 calculation algorithm Methods 0.000 claims description 11
- 230000001174 ascending effect Effects 0.000 claims description 7
- 230000009514 concussion Effects 0.000 claims description 6
- 238000000137 annealing Methods 0.000 claims 2
- 238000004064 recycling Methods 0.000 claims 1
- 238000007711 solidification Methods 0.000 abstract description 14
- 230000008023 solidification Effects 0.000 abstract description 14
- 238000010586 diagram Methods 0.000 description 12
- 238000004364 calculation method Methods 0.000 description 9
- 230000008859 change Effects 0.000 description 8
- 239000006185 dispersion Substances 0.000 description 8
- 230000008569 process Effects 0.000 description 6
- 230000003595 spectral effect Effects 0.000 description 6
- 238000005516 engineering process Methods 0.000 description 5
- 238000004611 spectroscopical analysis Methods 0.000 description 5
- 230000005012 migration Effects 0.000 description 4
- 238000013508 migration Methods 0.000 description 4
- 238000003860 storage Methods 0.000 description 4
- 238000013459 approach Methods 0.000 description 3
- 230000008901 benefit Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 239000003989 dielectric material Substances 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000009466 transformation 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
- 238000009825 accumulation Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000005266 casting Methods 0.000 description 1
- 230000006854 communication Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 230000004907 flux Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012913 prioritisation Methods 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 238000010845 search algorithm Methods 0.000 description 1
- 238000004513 sizing Methods 0.000 description 1
- 230000007480 spreading Effects 0.000 description 1
- 238000003892 spreading Methods 0.000 description 1
- 238000006467 substitution reaction 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
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (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)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明实施例提供一种全局优化的交错网格有限差分正演模拟方法和装置，其中，该方法包括：获得交错网格有限差分算子；基于最大化范数建立目标函数，求解目标函数得到用于优化交错网格有限差分算子的交错网格有限差分权系数，交错网格有限差分权系数是在优化后的交错网格有限差分算子的最大波数覆盖范围内且在优化后的交错网格有限差分算子的最大谱误差容限下计算得到的；采用求解得到的交错网格有限差分权系数优化空间交错网格有限差分算子并进行地震波正演模拟。该方案有利于更加灵活地控制优化后的交错网格有限差分算子的绝对谱误差，更有效地压制数值频散，提高地震模拟精度。
Description
技术领域
本发明涉及地震波数值模拟研究技术领域，特别涉及一种全局优化的交错网格有限差分正演模拟方法和装置。
背景技术
地震波正演模拟即模拟地震波在地球介质中的传播过程，并研究地震波的传播特征与地球介质参数的关系，达到与实际观测地震记录的最优逼近。由于实际地球介质的异常复杂，很难得到地震波波动方程的解析解。故地震波数值模拟方法是人们正确认识地震波的传播规律，验证各种假设地球模型正确性的主要方法技术；是地震勘探和地震学的重要基础，对实际地震资料的地质解释和储层预测等工作具有重要的理论指导作用和实际利用价值。由于广泛地被应用在地震勘探和天然地震领域中，地震数值模拟方法也得到迅速发展和有效应用。目前，地震波正演模拟的数值方法主要有：射线追踪法、积分方程法和波动方程法三大类。波动方程法是建立在地震波传播方程的理论基础上，对复杂介质中地震波传播具有广泛适应性，在地震波数值模拟中应用最广泛。
有限差分法是偏微分方程的主要数值解法之一，是最早用于地震波数值模拟中的数值模拟方法。波动方程有限差分法相比于其他数值模拟方法如伪谱法、有限元法，因为其存储的有效利用、计算量低且简单易于实现而更受人偏爱广泛应用在地震数值模拟、偏移和反演中。1968年Alterman和Karal最先尝试将有限差分法应用于层状介质弹性波传播的数值模拟中，产生了理论地震图。1972年，Boore将有限差分法用于非均匀介质地震波传播的模拟。随后，有限差分法被进一步用粘弹性、各向异性等复杂介质中和三维模型中地震波传播的模拟。
常规的有限差分法正演时，会产生严重的数值假频，称为数值频散，无法使波传播较大距离。而且时间域算法的误差会逐渐积累，导致计算精度不足，数值模拟结果分辨率降低。这是因为对模型进行离散成网格，并用差分算子逼近微分算子，从而使波动方程的系数发生变化即使相速度变成离散网格间隔的函数。因此当每个波长内采样点太少时，就会产生网格数值频散；如果增加单位波长内采样点，又会增加计算量和存储量，使正演成本增加。
如何消除这种数值频散，提高模拟精度是应用有限差分法正演的关键问题之一。为了达到这一目标，前人从各方面作了很多努力：在时间-空间域有限差分方面，从Alford等在1974年的二阶有限差分法到Dablain于1986年提出了高阶有限差分法算子，随着精度阶数增加，单位波长需要的采样点越来越少。一般情况下，时间方向采用二阶差分格式，空间方向采用高阶差分方式。从规则网格发展到交错网格(Virieux在1986首先提出；Ozdenvar&Mcmechan，1997；董良国，2000；KetilHokstad，2003)，交错网格方法提高了地震模拟的精确程度及计算稳定性，并消除部分假像。此外，Boris和Book1973年、1975提出用求解守恒方程的通量校正传输方法(FCT)，有效的压制了在粗网格上用差分算子引起的数值频散。
由于用高阶有限差分算子计算每个格点值时，参与计算的点数增加，使高阶有限差分的计算量大大增加。而且有人提出高阶有限差分算子存在饱和效应，且计算不稳定。极限有限差分法即伪谱法，不存在数值频散问题且每个波长只需要两个采样点，但计算量太大。
Zhou和Greenhalgh(1992)使用了广义的汉宁(Hanning)窗截断伪谱法算子得到了一种优化的有限差分算子；Igel et al.(1995)使用了高斯窗截断得到了一种优化的有限差分算子。优化的有限差分算子相对于常规有限差分算子其谱覆盖范围增大，即利用低阶的有限差分算子可以达到常规高阶算子的精度。Chu and stoffa2012利用二项式窗函数统一了有限差分法和伪谱法：用二项式窗函数对伪谱法算子进行截断，得到常规有限差分系数；并改进了二项式窗函数得到了优化的有限差分系数，有效提高了差分精度。wang等2014、2015年利用组合的窗函数直观可视化地控制其主瓣和旁瓣的形状，进而调节有限差分算子逼近微分算子的精度。Sun在2015年通过对一阶导数的有限宽频谱进行傅里叶(Fourier)反变换，得出了与chu在2012文中一致的交错网格差分算子形式，并用正规hanning窗函数截断后得到优化的差分系数。
Holberg于1987用最优化方法最小化群速度的相对峰值误差，而Etgen于2007提出用相速度构建优化问题的目标函数，由此均得到了优化的差分系数。Liu等2009提出同时考虑时间导数和空间导数，基于时空域频散关系，用泰勒展开法求解地震波方程的有限差分算子；并在2013、2014用最小二乘理论确定了的时空域优化的空间有限差分算子。Zhang等于2013年提出了在给定目前最小的误差限下(万分之一)，用模拟退火算法直接对近似空间二阶导数的有限差分算子进行优化，使其覆盖更大的波数范围。Yang、Yan等于2014、2015、2016年分别用最小二乘法、样点逼近法及雷米兹(Remez)算法确定了优化的空间交错网格有限差分算子。Liang等在2013、2014用线性化方法确定了时空域优化的空间有限差分算子。辛维等2015分别用线性化算子和模拟退火法确定了给定波数上限下优化的空间交错网格有限差分算子，并进行了对比。以上这些方法，都一定程度上减小了有限差分模拟的数值频散，并提高了数值模拟的计算效率。
综上所述，对有限差分算子进行进一步的优化，可分为窗函数优化算法和直接优化算法。窗函数优化方法是选择合适的窗函数去截断伪普法算子，从而得到优化的有限差分算子。窗函数优化法的一个关键问题就是窗函数的选择，其直接影响有限差分算子精度。该方法虽然优化过程直观，但是无法直接控制有限差分算子的精度，且受窗函数限制也无法达到很高的精度。直接优化算法是则是利用最小二乘法、Remez算法等最优化方法直接求使差分算子在波数域的误差曲线尽可能满足需求的差分系数。该类方法使有限差分算子的设计简化成一优化问题，给定谱覆盖范围，用优化方法搜索使误差尽可能小的有限差分系数。
针对现有技术的实现方案，Yang在2014年提出的基于最小二乘法理论的交错网格有限差分算子优化方案中，是用最小二乘法对给定波数范围内的误差的积分求极小，从而确定一种优化的差分系数。Yang在2015年提出的线性化方法求解优化交错网格有限差分算子方案中，也需要给定波数覆盖范围，并进行等间隔采样；带入误差函数后形成一线性方程组，然后用高斯消元法求解得唯一解。辛维于2015年应用最大最小原则作为模拟退火法优化算法求解交错网格有限差分算子的数值频散误差判定标准，得到了特定波数覆盖范围下优化的差分系数。
然而，上述两种方案以及Yang在2016年提出的Remez算法均存在以下缺陷：1)均需要预先设定优化差分算子所能覆盖的最大波数范围；2)均不能直接控制频散误差，需要先画出相对误差曲线才能进行误差控制；3)仅对相对误差进行了分析，以相对误差大小判断差分算子精度不够精确。因此，灵活性差、应用局限。
发明内容
本发明实施例提供了一种全局优化的交错网格有限差分正演模拟方法，以解决现有技术中优化交错网格有限差分算子的方案中存在的灵活性差、应用局限的技术问题。该方法包括：获得交错网格有限差分算子；基于最大化范数建立目标函数，求解所述目标函数得到用于优化所述交错网格有限差分算子的交错网格有限差分权系数，采用求解得到的交错网格有限差分权系数优化所述交错网格有限差分算子，得到优化的交错网格有限差分算子，其中，所述交错网格有限差分权系数是在优化的交错网格有限差分算子的最大波数覆盖范围内且在优化的交错网格有限差分算子的最大谱误差容限下计算得到的；利用所述优化的交错网格有限差分算子计算空间导数来进行地震波正演模拟。
在一个实施例中，所述目标函数为：
其中，E(c)为所述目标函数的极值；是优化的交错网格有限差分算子的最大波数覆盖范围；cn是第n个所述交错网格有限差分权系数；kx是波数覆盖范围；T是优化的交错网格有限差分算子的最大谱误差容限；N是得到所述交错网格有限差分算子所用的网格点数的长度；Δx是采样间隔。
在一个实施例中，求解所述目标函数，包括：采用模拟退火算法求解所述目标函数。
在一个实施例中，采用模拟退火算法求解所述目标函数，包括：设定所述优化的交错网格有限差分算子的最大谱误差容限为预设谱误差容限数值；通过循环以下步骤求解所述目标函数得到优化所述交错网格有限差分算子的相应交错网格有限差分权系数：将波数在0到尼奎斯特波数的范围内离散采样，按照波数采样点由小到大的顺序，先将最小的一个波数采样点作为预设波数范围的上限，该预设波数范围的下限为0，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数，在该得到的交错网格有限差分权系数使得所述目标函数的极值小于所述预设谱误差容限数值时，将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数，在得到的交错网格有限差分权系数使得所述目标函数的极值大于所述预设谱误差容限数值时，终止计算，将上一个波数采样点采用模拟退火算法计算得到的交错网格有限差分权系数输出，作为所述优化的交错网格有限差分算子的交错网格有限差分权系数；在计算得到的交错网格有限差分权系数使得所述目标函数的极值小于所述预设谱误差容限数值时，继续再将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数。
在一个实施例中，在求解所述目标函数得到优化所述交错网格有限差分算子的每个交错网格有限差分权系数之前，还包括：设定模拟退火算法的自变量维度和各维度自变量的取值范围，其中，自变量维度表示所述交错网格有限差分权系数的个数，各维度自变量的取值范围表示各个所述交错网格有限差分权系数的振幅绝对值范围，将自变量维度设定为得到所述交错网格有限差分算子所用的网格点数的长度的二分之一；将各维度自变量的取值范围设定为区间【0，2】；设定各个所述交错网格有限差分权系数的振幅围绕中心位置呈衰减的震荡。
本发明实施例还提供了一种全局优化的交错网格有限差分正演模拟装置，以解决现有技术中优化交错网格有限差分算子的方案中存在的灵活性差、应用局限的技术问题。该装置包括：有限差分算子获取模块，用于获得交错网格有限差分算子；目标函数求解模块，用于基于最大化范数建立目标函数，求解所述目标函数得到用于优化所述交错网格有限差分算子的交错网格有限差分权系数，采用求解得到的交错网格有限差分权系数优化所述交错网格有限差分算子，得到优化的交错网格有限差分算子，其中，所述交错网格有限差分权系数是在优化的交错网格有限差分算子的最大波数覆盖范围内且在优化的交错网格有限差分算子的最大谱误差容限下计算得到的；地震波正演模块，用于利用所述优化的交错网格有限差分算子计算空间导数来进行地震波正演模拟。
在一个实施例中，所述目标函数为：
其中，E(c)为所述目标函数的极值；是优化的交错网格有限差分算子的最大波数覆盖范围；cn是第n个所述交错网格有限差分权系数；kx是波数覆盖范围；T是优化的交错网格有限差分算子的最大谱误差容限；N是得到所述交错网格有限差分算子所用的网格点数的长度；Δx是采样间隔。
在一个实施例中，所述目标函数求解模块具体用于采用模拟退火算法求解所述目标函数。
在一个实施例中，所述目标函数求解模块，包括：误差设定单元，用于设定所述优化的交错网格有限差分算子的最大谱误差容限为预设谱误差容限数值；目标函数求解单元，用于通过循环以下步骤求解所述目标函数得到优化所述交错网格有限差分算子的相应交错网格有限差分权系数：将波数在0到尼奎斯特波数的范围内离散采样，按照波数采样点由小到大的顺序，先将最小的一个波数采样点作为预设波数范围的上限，该预设波数范围的下限为0，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数，在该得到的交错网格有限差分权系数使得所述目标函数的极值小于所述预设谱误差容限数值时，将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数，在得到的交错网格有限差分权系数使得所述目标函数的极值大于所述预设谱误差容限数值时，终止计算，将上一个波数采样点采用模拟退火算法计算得到的交错网格有限差分权系数输出，作为所述优化的交错网格有限差分算子的交错网格有限差分权系数；在得到的交错网格有限差分权系数使得所述目标函数的极值小于所述预设谱误差容限数值时，继续再将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数。
在一个实施例中，所述目标函数求解模块，还包括：参数设定单元，用于在求解所述目标函数得到优化所述交错网格有限差分算子的每个间交错网格有限差分权系数之前，设定模拟退火算法的自变量维度和各维度自变量的取值范围，其中，自变量维度表示所述交错网格有限差分权系数的个数，各维度自变量的取值范围表示各个所述交错网格有限差分权系数的振幅绝对值范围，将自变量维度设定为得到所述交错网格有限差分算子所用的网格点数的长度的二分之一；将各维度自变量的取值范围设定为区间【0，2】；设定各个所述交错网格有限差分权系数的振幅围绕中心位置呈衰减的震荡。
在本发明实施例中，通过基于最大化范数建立目标函数，并求解目标函数得到用于优化所述交错网格有限差分算子的交错网格有限差分权系数，采用求解得到的交错网格有限差分权系数优化交错网格有限差分算子，得到优化的交错网格有限差分算子，进而利用优化的交错网格有限差分算子计算空间导数来进行地震波正演模拟。由于目标函数是基于最大化范数建立的，相对于现有技术中最小平方法建立的目标函数，使得优化交错网格有限差分算子的过程更加灵活；同时，由于通过求解目标函数得到的交错网格有限差分权系数是在优化的交错网格有限差分算子的最大波数覆盖范围内且在优化的交错网格有限差分算子的最大谱误差容限下计算得到的，有利于更加灵活地控制优化的交错网格有限差分算子的绝对谱误差，使得优化的交错网格有限差分算子在较紧的谱误差容限内还能够覆盖尽可能大的波数范围，从而更有效地压制数值频散，提高地震模拟精度，有利于为逆时偏移、全波形反演提供高精度、高效的正演模拟方案，降低应用局限的问题。
附图说明
此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的限定。在附图中：
图1是本发明实施例提供的一种全局优化的交错网格有限差分正演模拟方法的流程图；
图2是本发明实施例提供的一种不同谱误差容限T下的8阶交错网格有限差分权系数绝对误差对比示意图；
图3是本发明实施例提供的一种图2的放大示意图；
图4是本发明实施例提供的一种采用模拟退火算法求解目标函数的流程图；
图5(a)是本发明实施例提供的一种谱误差容限T＝0.0001时的高阶优化的交错网格差分算子绝对误差对比示意图；
图5(b)是本发明实施例提供的一种图5(a)的放大示意图；
图6(a)是本发明实施例提供的一种700ms时刻8阶常规交错网格有限差分算子模拟X分量的示意图；
图6(b)是本发明实施例提供的一种700ms时刻8阶优化的交错网格有限差分算子模拟X分量的示意图；
图6(c)是本发明实施例提供的一种700ms时刻8阶常规交错网格有限差分算子模拟Z分量的示意图；
图6(d)是本发明实施例提供的一种700ms时刻8阶优化的交错网格有限差分算子模拟Z分量的示意图；
图7是本发明实施例提供的一种纵波速度模型示意图；
图8(a)是本发明实施例提供的一种marmousi2模型弹性波正演的垂直分量波场记录的示意图；
图8(b)是本发明实施例提供的一种图8(a)中区域1的放大示意图；
图8(c)是本发明实施例提供的一种图8(a)中区域2的放大示意图；
图9是本发明实施例提供的一种全局优化的交错网格有限差分正演模拟装置的结构框图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图，对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为对本发明的限定。
从现有技术中，我们发现用交错网格有限差分数值进行数值模拟时，由于网格离散导致的数值频散仍不可避免。关键是在保证精度的情况下，如何增大谱覆盖范围，以减少单位波长内采样点数，从而减少存储量和计算量；在单位波长采样点固定时，则能更好地压制高波数数值频散以适应高频地震模拟。基于泰勒展开的常规交错有限差分算子所覆盖的谱范围最小，但是波数较小时能最精确的逼近微分算子。基于最小二乘法优化的差分算子考虑奈奎斯特(Nyquist)波数内所有波数下的误差积分极小，结果所覆盖的谱范围最大；但是在波数较小时与微分算子的谱误差较大，无法有效压制数值频散。因此，压制交错网格有限差分算子的关键是在误差范围和谱覆盖范围之间找到一种平衡，本申请发明人提出上述全局优化的交错网格有限差分正演模拟方法，利用最大化范数建立的目标函数，则可以灵活地控制谱误差容限，使给定波数范围内所有样点的谱误差均在该谱误差容限内。因此，首先找到合适的优化算法求解该目标函数，然后选择误差容限，使得到的优化的交错网格有限差分算子在较紧的误差限内还能够覆盖尽可能大的波数范围。从而更有效地压制数值频散，为逆时偏移、全波形反演提供高精度、高效的正演模拟方案。
在本发明实施例中，提供了一种全局优化的交错网格有限差分正演模拟方法，如图1所示，该方法包括：
步骤101：获得交错网格有限差分算子；
步骤102：基于最大化范数建立目标函数，求解所述目标函数得到用于优化所述交错网格有限差分算子的交错网格有限差分权系数，采用求解得到的交错网格有限差分权系数优化所述交错网格有限差分算子，得到优化的交错网格有限差分算子，其中，所述交错网格有限差分权系数是在优化的交错网格有限差分算子的最大波数覆盖范围内且在优化的交错网格有限差分算子的最大谱误差容限下计算得到的；
步骤103：利用所述优化的交错网格有限差分算子计算空间导数来进行地震波正演模拟。
由图1所示的流程可知，在本发明实施例中，通过基于最大化范数建立目标函数，并求解目标函数得到用于优化所述交错网格有限差分算子的交错网格有限差分权系数，采用求解得到的交错网格有限差分权系数优化交错网格有限差分算子，得到优化的交错网格有限差分算子，进而利用优化的交错网格有限差分算子计算空间导数来进行地震波正演模拟。由于目标函数是基于最大化范数建立的，相对于现有技术中最小平方法建立的目标函数，使得优化交错网格有限差分算子的过程更加灵活；同时，由于通过求解目标函数得到的交错网格有限差分权系数是在优化的交错网格有限差分算子的最大波数覆盖范围内且在优化的交错网格有限差分算子的最大谱误差容限下计算得到的，有利于更加灵活地控制优化的交错网格有限差分算子的绝对谱误差，使得优化的交错网格有限差分算子在较紧的谱误差容限内还能够覆盖尽可能大的波数范围，从而更有效地压制数值频散，提高地震模拟精度，有利于为逆时偏移、全波形反演提供高精度、高效的正演模拟方案，降低应用局限的问题。
具体实施时，在进行地震波正演模拟的过程中，利用优化的交错网格有限差分算子近似计算空间导数，时间导数则可以用常规二阶交错网格有限差分算子近似计算。
具体实施时，本申请对获得上述交错网格有限差分算子(该交错网格有限差分算子可以为高阶交错网格有限差分算子)的方式不做限定，用sinc函数插值理论可以推导出高阶交错网格有限差分算子，当然用连续傅里叶反变换也可以推导出的一致的结果。例如，我们用sinc函数插值理论通过以下步骤获得高阶交错网格有限差分算子：
根据离散信号的采样理论，一个带限的连续信号f(x)可以被以一个均匀采样的信号fn通过sinc函数插值重建：
其中，Δx为采样间隔，为截止波数。
我们为得到函数f(x)的一阶导数，可以对公式(1)左右两边分别微分，即：
对于一阶导数，交错网格有限差分算子相对于常规有限差分算子具有更高的精度和更佳的稳定性，通常在波场模拟中获得更好的结果。因此，为推导出一阶导数的高阶交错网格有限差分算子，我们选择用公式(2)在处来估计一阶导数如下：
因此，用一个长度为N点的窗函数(其中N为偶数)，去截断公式(3)，可得到N点空间交错网格有限差分算子：
对于常规交错网格有限差分算子，为二项式窗函数(Chu和Stoffa，2012)，表示如下：
式中，l≤m；
是二项式系数公式，l≤m。也就是说，常规交错网格有限差分算子权系数对应于微分的sinc插值系数乘以加权的的二项式窗函数。对于一阶导数的交错网格有限差分权系数cn可以用加权的二项式窗函数找到，如下：
由公式(8)注意到，权系数因此公式(7)可以写成：
至此，我们基于sinc函数插值理论推导出了N点常规交错网格有限差分算子。
对公式(9)进行Fourier变换可得：
式中，kx为波数覆盖范围。此外，将公式(10)右边的正弦函数用泰勒级数展开，可建立N阶线性方程组，通过求解该方程组也能推导出交错网格常规有限差分算子。正因为常规交错网格有限差分算子无论是高阶还是低阶算子均是围绕零波数进行泰勒展开逼近，使得其在小波数范围内模拟精度很高，而在大波数范围内由于误差较大而导致宽频地震模拟中严重的空间频散。
具体实施时，为适应宽频带地震模拟需要，我们必须拓宽交错网格有限差分算子的波数覆盖范围，提高高频地震模拟精度，因此，在本实施例中，利用最大化范数控制一定波数范围内的最大绝对谱误差，建立上述目标函数如下：
其中，E(c)为所述目标函数的极值；是优化的交错网格有限差分算子的最大波数覆盖范围；cn是第n个所述交错网格有限差分权系数；kx是波数覆盖范围；T是优化的交错网格有限差分算子的最大谱误差容限；N是得到所述交错网格有限差分算子所用的网格点数的长度；Δx是采样间隔。使得我们可以得到谱覆盖范围更大的优化的交错网格有限差分权系数。很明显，基于最大化范数建立的目标函数比最小平方法建立的目标函数更加灵活。
具体实施时，在建立目标函数之后，可知谱误差容限T决定了交错网格有限差分算子能否寻优成功以及其最终逼近精度。图2为谱误差容限T分别设为0.00005、0.0001、0.0002、0.0004和0.0005时，相应得到的8阶交错网格有限差分权系数的绝对误差谱曲线，并与常规交错网格差分系数(图2中实线曲线表示常规交错网格差分系数，实线曲线上的4、6、8、10等数字表示常规交错网格差分系数的阶数)进行对比。从图2中可以发现随着谱误差容限T逐渐变松，交错网格有限差分权系数的谱覆盖范围逐渐增大，T＝0.0005时，交错网格有限差分权系数最大且超过了16阶常规交错网格有限差分系数；从图3中则发现，随着谱误差容限T逐渐变松，交错网格有限差分权系数的谱误差在小波数时波动逐渐变大。因此，我们可以得出结论：设定谱误差容限T若太小，则优化的交错网格有限差分算子覆盖的谱范围过小使宽频带弹性波模拟精度偏低；若设定谱误差容限T太大，则优化的交错网格有限差分算子覆盖谱范围较大但是误差积累影响比较严重，无法适应长时间步、深层弹性波模拟。
模拟退火法(SA)于1953年Metropolis提出属于一种随机算法，并由Kirkpatrick在1983年将其应用到组合优化问题。模拟退火算法采用随机采样方式寻找新解，并以随温度下降而逐渐变小的概率接受较差的解。因此，模拟退火算法是局部搜索算法的扩展，具有能够跳出局部最优陷阱的优点，较快的收敛到全局最优解。用模拟退火法寻优时初始温度和终止温度的选择非常重要，直接决定搜索多长，以及能否成功搜索到满足误差上限的解。因此，针对如此复杂的目标函数，线性搜索方法和最小二乘法无法灵活求解。本申请提出采用模拟退火法对上述目标函数进行求解，直接搜索最优的交错网格有限差分权系数。
具体实施时，采用模拟退火算法求解所述目标函数，包括：设定所述优化的交错网格有限差分算子的最大谱误差容限为预设谱误差容限数值；通过循环以下步骤求解所述目标函数得到优化所述交错网格有限差分算子的相应交错网格有限差分权系数：将波数在0到尼奎斯特波数的范围内离散采样，按照波数采样点由小到大的顺序，先将最小的一个波数采样点作为预设波数范围的上限，该预设波数范围的下限为0，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数，在该得到的交错网格有限差分权系数使得所述目标函数的极值小于所述预设谱误差容限数值时，将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数，在得到的交错网格有限差分权系数使得所述目标函数的极值大于所述预设谱误差容限数值时，终止计算，将上一个波数采样点采用模拟退火算法计算得到的交错网格有限差分权系数输出，作为所述优化的交错网格有限差分算子的交错网格有限差分权系数；在计算得到的交错网格有限差分权系数使得所述目标函数的极值小于所述预设谱误差容限数值时，继续再将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数。
具体的，本申请采用模拟退火算法求解上述目标函数，并不直接设定优化的交错网格有限差分算子所能覆盖的最大波数范围，而是在0到尼奎斯特波数的范围内离散采样后，在波数采样点为最大波数上限的预设波数范围中，搜索优化的交错网格有限差分算子在满足谱误差容限数值上限的情况下所能覆盖的最大波数范围。
具体实施时，为了减少优化代价，在本实施例中，在求解所述目标函数得到优化所述交错网格有限差分算子的每个交错网格有限差分权系数之前，还包括：设定模拟退火算法的自变量维度和各维度自变量的取值范围，其中，自变量维度表示所述交错网格有限差分权系数的个数，各维度自变量的取值范围表示各个所述交错网格有限差分权系数的振幅绝对值范围，将自变量维度设定为得到所述交错网格有限差分算子所用的网格点数(即窗函数)的长度的二分之一；将各维度自变量的取值范围设定为区间【0，2】；设定各个所述交错网格有限差分权系数的振幅围绕中心位置呈衰减的震荡。
具体的，本申请发明人根据sinc插值理论发现：(1)优化的N点交错网格有限差分算子中有N/2个交错网格有限差分权系数，比如N＝16，则有8个系数要优化；(2)各个交错网格有限差分权系数的振幅绝对值均在[0，2]区间内；(3)交错网格有限差分权系数的振幅应该是围绕中心位置0呈衰减的震荡，即|cn|＞|cn+1|，且|cn||cn+1|＜0，for n＝1,2,...,N/2。这样我们可以通过只确定c1到cn而优化整个交错网格有限差分算子。在设定模拟退火算法的关键参数时，可以根据上述(1)和(2)模拟退火算法的自变量维度M和各维度自变量取值范围，再用(3)原则对解进行限制大大减少了搜索范围，使模拟退火法更高效。事实上通过观察传统差分算子，发现原始系数也满足这三点约束。这就表明这三个规则对通常的差分算子都是合理、适用的。
具体的，采用模拟退火算法求解上述目标函数并进行地震波传播模拟的过程如图4所示：
(1)、优化的交错网格有限差分算子的最大谱误差容限T为预设谱误差容限；
(2)、将波数在范围内均匀离散采样，按照波数采样点由小到大的顺序，K表示波数采样点的序数，是以一个波数采样点作为上限的预设波数范围，从K＝1开始对每个波数采样点执行第3-4步，搜索满足T的最优交错网格有限差分权系数。
(3)、用模拟退火算法进行全局搜索，其中要设定关键参数有：马尔可夫(markov)链长度、起始温度、终止温度和步长因子；此外，还有自变量维度M和各维度自变量取值范围，即2M阶交错网格有限差分算子要优化的交错网格有限差分系数个数及各交错网格有限差分权系数的取值范围。
(4)、若搜索到有交错网格有限差分权系数使目标函数的极值小于预设谱误差容限，则回到第二步，且k＝k+1，即将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解目标函数，进行下一个波数采样点的全局寻优；若使目标函数的极值大于预设谱误差容限，则终止，取上一个波数采样点的交错网格有限差分系数作为最优解输出。
(5)、利用2M阶优化的交错网格有限差分算子计算空间导数来进行地震波传播的高精度模拟。
优化的交错网格有限差分算子更佳地逼近一阶空间微分算子，1)可以用于各种介质中的高精度声波模拟和逆时偏移成像，包括各向同性介质和各向异性介质；2)可以用于各种介质中的高精度弹性波模拟和逆时偏移成像；3)对于基于CPU/GPU协同并行的地震波模拟装置，仅需要用优化的交错网格差分权系数代替原来的差分系数即可得到数值频散更小的模拟结果。其中以弹性波的模拟为应用实例，给出二维情况下、各向同性介质中的一阶交错网格速度-应力方程如下：
式中，λ、μ是两个拉梅常数，
具体实施时，在应用上述全局优化的交错网格有限差分正演模拟方法时，我们设定较紧的谱误差容限T为0.0001，即以保证波在传播很长距离后积累误差仍然很小可以接受。然后采用模拟退火算法，求解交错网格有限差分权系数，较快地收敛到全局最优解，由此得到优化的交错网格有限差分算子：(1)使谱绝对误差控制在万分之一内，这是目前所达到的最小误差；(2)尽可能地覆盖更大的谱范围，这也是目前此谱误差容限下所达到的最大谱覆盖范围。从而在利用优化的交错网格有限差分算子进行波场模拟过程中，可以有效地降低空间离散引入的频散误差，提高模拟结果的精度。因此，我们称本专利优化的交错网格有限差分算子相比于相同阶数下的常规网格有限差分算子，具有更高的精度。从图5(a)(图5(a)中实线曲线表示常规交错网格差分系数，实线曲线上的4、6、8、10等数字表示常规交错网格差分系数的阶数)中可以看出，优化的8阶交错网格有限差分算子达到常规12阶交错网格差分算子的精度；优化的12阶交错网格有限差分算子则超过了常规24阶交错网格差分算子的精度。从图5(b)中可以发现，优化的8阶交错网格有限差分算子在精度误差0.0001下的谱覆盖范围远超过16阶常规交错网格差分算子；而优化的10、12阶交错网格有限差分算子在精度误差0.0001下的谱覆盖范围远超过常规28阶交错网格差分算子。因此，用低阶优化的交错网格有限差分算子代替高阶常规交错网格有限差分算子，由于算子长度变短而减少了计算量和计算时间；另外，基于GPU并行计算平台，差分算子长度越短则对应全局内存访问量越少，减少延迟时间。即本专利优化的交错网格有限差分算子在保证了模拟精度的条件下，还可以有效地提高地震波模拟的效率。
以下结合具体示例来说明上述全局优化的交错网格有限差分正演模拟方法的有效性。
实例一、将优化的8阶交错网格有限差分算子近似计算空间一阶导数，而时间一阶导数仍用常规二阶交错网格差分算子计算，求解公式(12)，模拟弹性波在均匀各向同性介质中的传播，以验证其有效性。
速度模型网格大小为500×500，横向和纵向网格间隔均为10m，纵波速度为2500m/s，横波速度为1700m/s，密度为常数1g/cm3；采用主频为33Hz的雷克(Ricker)子波点源加载在垂向应力上，在模型中心激发，时间采样间隔为1ms。如图6(a)、6(b)、6(c)以及6(d)所示，可以发现在算子长度相同时(相同阶数)，优化的交错网格有限差分算子模拟纵波和横波的频散明显要比常规交错网格差分算子模拟的少。从而，证明了本专利方案压制数值频散、提高地震波模拟精度方面的有效性。
实例二、用复杂的marmousi2模型验证优化的8阶交错网格有限差分算子的压制频散性能。
Marmousi2模型大小为961*561，其纵波速度模型如图7所示，横波速度用给出，密度为常数1g/cm3。采用主频为38Hz的Ricker子波，在靠近地表处(2.9976，0.01249)Km激发；横向和纵向网格间隔均为6.245m，时间采样间隔为0.5ms。图8(a)为优化的8阶交错网格有限差分方法和常规方法的时间总长度为6s的垂直分量波场记录，其中，左侧是8阶常规交错网格有限差分模拟；右侧为8阶优化交错网格有限差分模拟；图8(b)、8(c)分别是图8(a)中区域1和2的放大图，可以看出优化的交错网格有限差分方法有效地压制了marmousi2中由横波低速引起的数值频散，提高了宽频弹性波模拟的精度。
基于同一发明构思，本发明实施例中还提供了一种全局优化的交错网格有限差分正演模拟装置，如下面的实施例所述。由于全局优化的交错网格有限差分正演模拟装置解决问题的原理与全局优化的交错网格有限差分正演模拟方法相似，因此全局优化的交错网格有限差分正演模拟装置的实施可以参见全局优化的交错网格有限差分正演模拟方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。
图9是本发明实施例的全局优化的交错网格有限差分正演模拟装置的一种结构框图，如图9所示，该装置包括：
有限差分算子获取模块901，用于获得交错网格有限差分算子；
目标函数求解模块902，用于基于最大化范数建立目标函数，求解所述目标函数得到用于优化所述交错网格有限差分算子的交错网格有限差分权系数，采用求解得到的交错网格有限差分权系数优化所述交错网格有限差分算子，得到优化的交错网格有限差分算子，其中，所述交错网格有限差分权系数是在优化的交错网格有限差分算子的最大波数覆盖范围内且在优化的交错网格有限差分算子的最大谱误差容限下计算得到的；
地震波正演模块903，用于利用所述优化的交错网格有限差分算子计算空间导数来进行地震波正演模拟。
在一个实施例中，所述目标函数为：
其中，E(c)为所述目标函数的极值；是优化的交错网格有限差分算子的最大波数覆盖范围；cn是第n个所述交错网格有限差分权系数；kx是波数覆盖范围；T是优化的交错网格有限差分算子的最大谱误差容限；N是得到所述交错网格有限差分算子所用的网格点数的长度；Δx是采样间隔。
在一个实施例中，所述目标函数求解模块具体用于采用模拟退火算法求解所述目标函数。
在一个实施例中，所述目标函数求解模块，包括：误差设定单元，用于设定所述优化的交错网格有限差分算子的最大谱误差容限为预设谱误差容限数值；目标函数求解单元，用于通过循环以下步骤求解所述目标函数得到优化所述交错网格有限差分算子的相应交错网格有限差分权系数：将波数在0到尼奎斯特波数的范围内离散采样，按照波数采样点由小到大的顺序，先将最小的一个波数采样点作为预设波数范围的上限，该预设波数范围的下限为0，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数，在该得到的交错网格有限差分权系数使得所述目标函数的极值小于所述预设谱误差容限数值时，将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数，在得到的交错网格有限差分权系数使得所述目标函数的极值大于所述预设谱误差容限数值时，终止计算，将上一个波数采样点采用模拟退火算法计算得到的交错网格有限差分权系数输出，作为所述优化的交错网格有限差分算子的交错网格有限差分权系数；在得到的交错网格有限差分权系数使得所述目标函数的极值小于所述预设谱误差容限数值时，继续再将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数。
在一个实施例中，所述目标函数求解模块，还包括：参数设定单元，用于在求解所述目标函数得到优化所述交错网格有限差分算子的每个交错网格有限差分权系数之前，设定模拟退火算法的自变量维度和各维度自变量的取值范围，其中，自变量维度表示所述交错网格有限差分权系数的个数，各维度自变量的取值范围表示各个所述交错网格有限差分权系数的振幅绝对值范围，将自变量维度设定为得到所述交错网格有限差分算子所用的网格点数的长度的二分之一；将各维度自变量的取值范围设定为区间【0，2】；设定各个所述交错网格有限差分权系数的振幅围绕中心位置呈衰减的震荡。
在本发明实施例中，通过基于最大化范数建立目标函数，并求解目标函数得到用于优化所述交错网格有限差分算子的交错网格有限差分权系数，采用求解得到的交错网格有限差分权系数优化交错网格有限差分算子，得到优化的交错网格有限差分算子，进而利用优化的交错网格有限差分算子计算空间导数来进行地震波正演模拟。由于目标函数是基于最大化范数建立的，相对于现有技术中最小平方法建立的目标函数，使得优化交错网格有限差分算子的过程更加灵活；同时，由于通过求解目标函数得到的交错网格有限差分权系数是在优化的交错网格有限差分算子的最大波数覆盖范围内且在优化的交错网格有限差分算子的最大谱误差容限下计算得到的，有利于更加灵活地控制优化的交错网格有限差分算子的绝对谱误差，使得优化的交错网格有限差分算子在较紧的谱误差容限内还能够覆盖尽可能大的波数范围，从而更有效地压制数值频散，提高地震模拟精度，有利于为逆时偏移、全波形反演提供高精度、高效的正演模拟方案，降低应用局限的问题。
显然，本领域的技术人员应该明白，上述的本发明实施例的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，并且在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明实施例不限制于任何特定的硬件和软件结合。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (8)
1.一种全局优化的交错网格有限差分正演模拟方法，其特征在于，包括：
获得交错网格有限差分算子；
基于最大化范数建立目标函数，求解所述目标函数得到用于优化所述交错网格有限差分算子的交错网格有限差分权系数，采用求解得到的交错网格有限差分权系数优化所述交错网格有限差分算子，得到优化的交错网格有限差分算子，其中，所述交错网格有限差分权系数是在优化的交错网格有限差分算子的最大波数覆盖范围内且在优化的交错网格有限差分算子的最大谱误差容限下计算得到的；
利用所述优化的交错网格有限差分算子计算空间导数来进行地震波正演模拟；
所述目标函数为：
其中，E(c)为所述目标函数的极值；是优化的交错网格有限差分算子的最大波数覆盖范围；cn是第n个所述交错网格有限差分权系数；kx是波数覆盖范围；T是优化的交错网格有限差分算子的最大谱误差容限；N是得到所述交错网格有限差分算子所用的网格点数的长度；Δx是采样间隔。
2.如权利要求1所述的全局优化的交错网格有限差分正演模拟方法，其特征在于，求解所述目标函数，包括：
采用模拟退火算法求解所述目标函数。
3.如权利要求2所述的全局优化的交错网格有限差分正演模拟方法，其特征在于，采用模拟退火算法求解所述目标函数，包括：
设定所述优化的交错网格有限差分算子的最大谱误差容限为预设谱误差容限数值；
通过循环以下步骤求解所述目标函数得到优化所述交错网格有限差分算子的相应交错网格有限差分权系数：
将波数在0到尼奎斯特波数的范围内离散采样，按照波数采样点由小到大的顺序，先将最小的一个波数采样点作为预设波数范围的上限，该预设波数范围的下限为0，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数，在该得到的交错网格有限差分权系数使得所述目标函数的极值小于所述预设谱误差容限数值时，将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数，在得到的交错网格有限差分权系数使得所述目标函数的极值大于所述预设谱误差容限数值时，终止计算，将上一个波数采样点采用模拟退火算法计算得到的交错网格有限差分权系数输出，作为所述优化的交错网格有限差分算子的交错网格有限差分权系数；在计算得到的交错网格有限差分权系数使得所述目标函数的极值小于所述预设谱误差容限数值时，继续再将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数。
4.如权利要求3所述的全局优化的交错网格有限差分正演模拟方法，其特征在于，在求解所述目标函数得到优化所述交错网格有限差分算子的每个交错网格有限差分权系数之前，还包括：
设定模拟退火算法的自变量维度和各维度自变量的取值范围，其中，自变量维度表示所述交错网格有限差分权系数的个数，各维度自变量的取值范围表示各个所述交错网格有限差分权系数的振幅绝对值范围，将自变量维度设定为得到所述交错网格有限差分算子所用的网格点数的长度的二分之一；将各维度自变量的取值范围设定为区间【0，2】；设定各个所述交错网格有限差分权系数的振幅围绕中心位置呈衰减的震荡。
5.一种全局优化的交错网格有限差分正演模拟装置，其特征在于，包括：
有限差分算子获取模块，用于获得交错网格有限差分算子；
目标函数求解模块，用于基于最大化范数建立目标函数，求解所述目标函数得到用于优化所述交错网格有限差分算子的交错网格有限差分权系数，采用求解得到的交错网格有限差分权系数优化所述交错网格有限差分算子，得到优化的交错网格有限差分算子，其中，所述交错网格有限差分权系数是在优化的交错网格有限差分算子的最大波数覆盖范围内且在优化的交错网格有限差分算子的最大谱误差容限下计算得到的；
地震波正演模块，用于利用所述优化的交错网格有限差分算子计算空间导数来进行地震波正演模拟；
所述目标函数为：
其中，E(c)为所述目标函数的极值；是优化的交错网格有限差分算子的最大波数覆盖范围；cn是第n个所述交错网格有限差分权系数；kx是波数覆盖范围；T是优化的交错网格有限差分算子的最大谱误差容限；N是得到所述交错网格有限差分算子所用的网格点数的长度；Δx是采样间隔。
6.如权利要求5所述的全局优化的交错网格有限差分正演模拟装置，其特征在于，所述目标函数求解模块具体用于采用模拟退火算法求解所述目标函数。
7.如权利要求6所述的全局优化的交错网格有限差分正演模拟装置，其特征在于，所述目标函数求解模块，包括：
误差设定单元，用于设定所述优化的交错网格有限差分算子的最大谱误差容限为预设谱误差容限数值；
目标函数求解单元，用于通过循环以下步骤求解所述目标函数得到优化所述交错网格有限差分算子的相应交错网格有限差分权系数：
将波数在0到尼奎斯特波数的范围内离散采样，按照波数采样点由小到大的顺序，先将最小的一个波数采样点作为预设波数范围的上限，该预设波数范围的下限为0，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数，在该得到的交错网格有限差分权系数使得所述目标函数的极值小于所述预设谱误差容限数值时，将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数，在得到的交错网格有限差分权系数使得所述目标函数的极值大于所述预设谱误差容限数值时，终止计算，将上一个波数采样点采用模拟退火算法计算得到的交错网格有限差分权系数输出，作为所述优化的交错网格有限差分算子的交错网格有限差分权系数；在得到的交错网格有限差分权系数使得所述目标函数的极值小于所述预设谱误差容限数值时，继续再将下一个波数采样点作为预设波数范围的上限，在该预设波数范围内采用模拟退火算法求解所述目标函数得到交错网格有限差分权系数。
8.如权利要求7所述的全局优化的交错网格有限差分正演模拟装置，其特征在于，所述目标函数求解模块，还包括：
参数设定单元，用于在求解所述目标函数得到优化所述交错网格有限差分算子的每个交错网格有限差分权系数之前，设定模拟退火算法的自变量维度和各维度自变量的取值范围，其中，自变量维度表示所述交错网格有限差分权系数的个数，各维度自变量的取值范围表示各个所述交错网格有限差分权系数的振幅绝对值范围，将自变量维度设定为得到所述交错网格有限差分算子所用的网格点数的长度的二分之一；将各维度自变量的取值范围设定为区间【0，2】；设定各个所述交错网格有限差分权系数的振幅围绕中心位置呈衰减的震荡。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710253444.XA CN106842306B (zh) | 2017-04-18 | 2017-04-18 | 一种全局优化的交错网格有限差分正演模拟方法和装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710253444.XA CN106842306B (zh) | 2017-04-18 | 2017-04-18 | 一种全局优化的交错网格有限差分正演模拟方法和装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106842306A CN106842306A (zh) | 2017-06-13 |
| CN106842306B true CN106842306B (zh) | 2019-03-01 |
Family
ID=59146950
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710253444.XA Expired - Fee Related CN106842306B (zh) | 2017-04-18 | 2017-04-18 | 一种全局优化的交错网格有限差分正演模拟方法和装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106842306B (zh) |
Families Citing this family (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107179549B (zh) * | 2017-07-11 | 2019-02-26 | 中海石油(中国)有限公司 | 一种时间域声波方程显式有限差分地震响应模拟方法 |
| CN107479092B (zh) * | 2017-08-17 | 2019-02-12 | 电子科技大学 | 一种基于方向导数的频率域高阶声波方程正演模拟方法 |
| CN109541686A (zh) * | 2017-09-22 | 2019-03-29 | 长江大学 | 具有普适性的频率域三维声波有限差分优化方法 |
| CN109116418B (zh) * | 2018-08-01 | 2020-07-10 | 中国石油天然气股份有限公司 | 三维波动方程混合网格有限差分数值模拟方法及装置 |
| CN109143340B (zh) * | 2018-08-20 | 2020-03-10 | 中国海洋石油集团有限公司 | 一种基于常q模型的粘弹介质地震波模拟方法及系统 |
| CN109490954B (zh) * | 2018-09-20 | 2019-12-20 | 中国科学院地质与地球物理研究所 | 波场正演模拟方法及装置 |
| CN109239776B (zh) * | 2018-10-16 | 2021-02-09 | 中国石油天然气股份有限公司 | 一种地震波传播正演模拟方法和装置 |
| CN109684760B (zh) * | 2018-12-29 | 2021-04-13 | 北京化工大学 | 基于随机搜索算法的弹性矢量波场数值模拟方法及系统 |
| CN109598093B (zh) * | 2018-12-29 | 2020-12-04 | 北京化工大学 | 基于拟合窗函数的地震矢量波场数值模拟方法及系统 |
| CN109711072B (zh) * | 2018-12-29 | 2023-04-07 | 北京化工大学 | 基于混合群智能算法的地震矢量波场数值模拟方法及系统 |
| CN110261896B (zh) * | 2019-04-26 | 2021-07-20 | 中国石油化工股份有限公司 | 稳定的各向异性ti介质正演模拟方法 |
| CN112379422A (zh) * | 2020-10-30 | 2021-02-19 | 中国石油天然气集团有限公司 | 垂变网格地震波场外推方法及装置 |
| US20240104097A1 (en) * | 2022-09-21 | 2024-03-28 | Ocient Holdings LLC | Caching precomputed binomial coefficient values for query execution |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103823239A (zh) * | 2013-10-13 | 2014-05-28 | 中国石油集团西北地质研究所 | 频率域优化混合交错网格有限差分正演模拟方法 |
| CN103630933B (zh) * | 2013-12-09 | 2017-01-18 | 中国石油天然气集团公司 | 基于非线性优化的时空域交错网格有限差分方法和装置 |
| CN103699798B (zh) * | 2013-12-25 | 2017-01-25 | 中国石油天然气股份有限公司 | 一种实现地震波场数值模拟方法 |
| CN104459773A (zh) * | 2014-08-08 | 2015-03-25 | 中国石油天然气集团公司 | 基于交错网格Lowrank分解的无条件稳定地震波场延拓方法 |
| CN105911584B (zh) * | 2015-09-25 | 2017-05-03 | 中国科学院地质与地球物理研究所 | 一种隐式交错网格有限差分弹性波数值模拟方法及装置 |
-
2017
- 2017-04-18 CN CN201710253444.XA patent/CN106842306B/zh not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| CN106842306A (zh) | 2017-06-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106842306B (zh) | 一种全局优化的交错网格有限差分正演模拟方法和装置 | |
| CN109923440B (zh) | 面波勘探方法及终端设备 | |
| Zhang et al. | Optimized explicit finite-difference schemes for spatial derivatives using maximum norm | |
| Druskin et al. | Near-optimal perfectly matched layers for indefinite Helmholtz problems | |
| Schwarz et al. | The use of FFT techniques in physical geodesy | |
| CN107153216B (zh) | 确定地震波场的坡印廷矢量的方法、装置以及计算机存储介质 | |
| He et al. | Determining the optimal coefficients of the explicit finite-difference scheme using the Remez exchange algorithm | |
| CN104133241B (zh) | 波场分离方法和装置 | |
| CN108196303B (zh) | 弹性波波场分离方法、装置、存储介质及设备 | |
| CN104459773A (zh) | 基于交错网格Lowrank分解的无条件稳定地震波场延拓方法 | |
| CN108919356B (zh) | 一种稳定的大沙漠区衰减补偿逆时偏移成像系统及方法 | |
| CN111766628A (zh) | 一种预条件的时间域弹性介质多参数全波形反演方法 | |
| Liu et al. | 3D acoustic wave modelling with time-space domain dispersion-relation-based finite-difference schemes and hybrid absorbing boundary conditions | |
| CN109541681A (zh) | 一种拖缆地震数据和少量obs数据联合的波形反演方法 | |
| CN107798156A (zh) | 一种频率域2.5维粘弹性波数值模拟方法及装置 | |
| Gadylshin et al. | Optimization of the training dataset for numerical dispersion mitigation neural network | |
| US20150134308A1 (en) | Method and device for acquiring optimization coefficient, and related method and device for simulating wave field | |
| WO2015155597A2 (en) | Attenuating pseudo s-waves in acoustic anisotropic wave propagation | |
| Xie et al. | Cosine-type weighted hybrid absorbing boundary based on the second-order Higdon boundary condition and its GPU implementation | |
| Jing et al. | A highly efficient time-space-domain optimized method with Lax-Wendroff type time discretization for the scalar wave equation | |
| Wang et al. | Time-varying multiscale spatial correlation: Simulation and application to wind loading of structures | |
| CN108363097A (zh) | 一种地震资料偏移成像方法 | |
| Fang et al. | Seismic wavefield modeling based on time-domain symplectic and Fourier finite-difference method | |
| Jing et al. | An optimized time-space-domain finite difference method with piecewise constant interpolation coefficients for scalar wave propagation | |
| Chen* et al. | Optimized time-space domain staggered-grid finite-difference methods based on new finite-difference stencils |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190301 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |