一种空间数值频散模拟和有限差分参数优化选取方法与流程

1.本发明涉及一种地震勘探有限差分正演数值模拟领域，具体地说它是一种有限差分数值模拟频散和参数优化选取方法。更具体地说它是一种空间数值频散模拟和有限差分参数优化选取方法。


背景技术：

2.有限差分方法具有实施方便、计算效率高等优点，被广泛应用于地震波场的数值模拟。数值频散问题一直是影响有限差分数值模拟精度的因素之一。有限差分数值模拟频散指，由于对微分方程进行差分离散数值计算，导致不同频率波场的传播速度不同。频率越大的波场，其数值模拟的相速度和群速度之间的误差也越大，频散现象也越严重。压制有限差分空间方向数值频散最直接的方法是缩小网格间距，提高网格划分精度，相对应也会使增加计算量。针对精细网格划分导致巨大计算代价问题和粗糙网格导致数值模拟频散问题，目前主要方法是优化高阶差分算子系数，旨在节约计算成本的同时可以高效率、高精度数值模拟出波场。例如：dablain(dablain m a.1896.the application of high-order differencing to the scalar wave equation.geophysics,51:54-66,doi:https://doi.org/10.1190/1.1442040)提出了基于泰勒级数展开方法的高阶精度有限差分数值模拟，holberg(holberg o.1987.computational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena.geophysical prospecting,35(6):629-655,doi:https://doi.org/10.1111/j.1365-2478.1987.tb00841.x.)提出了基于优化差分系数压制数值频散的思想，zhang(zhang j h,yao z x.2013.optimized finite-difference operator for broadband seismic wave modeling:geophysics,78(1):a13-a18,doi:10.1190/geo2012-0277.1)给出了相速度误差阈值为0.0001时，优化类算法求取出的高阶算子差分系数。liu(liu y.2013.globally optimal finite-difference schemes based on least squares.geophysics,78(4)：t113-t132,doi：10.1190/geo2012-0480.1.liu y.2020.acoustic and elastic finite-difference modeling by optimal variable-length spatial operators.geophysics，85(2)：t57-t70，doi：https://doi.org/10.1190/geo2019-0145.1)给出了基于最小二乘方法(least-squares method，ls)和雷米兹交换方法(remez exchange method，re)求取的高阶算子差分系数。
3.优化高阶算子差分系数是基于人为给定一个相速度或者群速度误差阈值，其差分系数可以一定程度上压制高频波场，但对于一定频带宽度波场的有限差分数值模拟时，无法为选取优化的有限差分参数提供指导。优化选取的有限差分参数可以保证在一定数值模拟精度的条件下，使有限差分数值模拟的效率最高，节省计算成本，促进生产效益。
4.因此，开发一种有限差分参数的优化选取方法很有必要。


技术实现要素：

5.本发明的目的是为了提供一种空间数值频散模拟和有限差分参数优化选取方法，为一种可以给出波场空间方向数值频散的解析表达式，并基于该解析表达式计算的空间频散误差，进行有限差分参数优化选取的方法。
6.为了实现上述目的，本发明的技术方案为：一种空间数值频散模拟和有限差分参数优化选取方法，其特征在于：包括如下步骤，
7.步骤一：给定震源波场；
8.步骤二：空间偏导数进行有限差分离散，得到映射波数，映射波数用于步骤三中；不同差分格式(例如：空间一阶、二偏导数的2n阶精度中心差分离散格式，空间一阶偏导数的2n阶精度交错网格差分离散格式等)对应着不同的波数映射关系及相对应的映射波数；
9.步骤三：对震源波场进行空间频散正变换(fsdt计算)，对没有频散的参考波场，进行空间频散正变换，得到有数值频散的波场，用于步骤四和步骤五；
10.步骤四：设定频散误差评定准则，用数值频散波场与无频散的震源波场之间误差作为评判频散程度的依据，不是真实波场，选取合适误差阈值，用于步骤五中筛选有限差分参数组合；
11.步骤五：依据给定合适误差阈值，筛选有限差分参数组合；
12.步骤六：设定目标准则，优化选取有限差分参数。
13.在上述技术方案中，在步骤一中，给定震源波场的具体方法为：选取理论震源子波或者根据真实地震记录提取出震源子波。
14.步骤一选取的震源子波为具有解析表达式的理论子波(例如高斯函数，雷克子波的一阶导数子波等)或者野外地震资料提取出来的子波。
15.在上述技术方案中，在步骤二中，得到映射波数的具体方法为：
16.步骤2.1：对空间偏导数进行有限差分离散(如，空间一阶、二偏导数的2n阶精度中心差分离散格式以及空间一阶偏导数的2n阶精度交错网格差分离散格式，其他差分格式参考类似推导)，获取差分离散表达式；
17.步骤2.2：对步骤2.1中的差分离散表达式进行傅里叶变换及简化，得到波数k的映射关系及映射波数k
′
。
18.在上述技术方案中，空间偏导数包括空间一阶偏导数、空间二阶偏导数。
19.在上述技术方案中，当空间偏导数选用空间一阶偏导数时，在步骤2.1中，获取空间一阶偏导数的2n阶精度中心差分离散表示为：
[0020][0021]
式(1)中：p(x,t)为关于空间坐标x以及时间变量t的波场变量；p
i l
表示t时刻位于空间离散位置(i l)δx处的波场；p
i-l
表示t时刻位于空间离散位置(i-l)δx处的波场；δx为空间离散步长；c
l
为高阶差分算子系数；
[0022]
空间一阶偏导数的2n阶精度中心差分离散波数映射关系及映射波数k
′
：
[0023][0024]
式(2)中：k为波数；
[0025]
空间一阶偏导数的2n阶精度交错网格差分离散可以表示为：
[0026][0027]
式(3)中：p(x,t)为关于空间坐标x以及时间变量t的波场变量；pi表示t时刻位于空间离散位置iδx处的波场；p
i l-1/2
表示t时刻位于空间离散位置(i l-1/2)δx处的波场；p
i-l 1/2
表示t时刻位于空间离散位置(i-l 1/2)δx处的波场；δx为空间离散步长；c
l
为高阶差分算子系数；
[0028]
空间一阶偏导数的2n阶精度交错网格差分离散波数映射关系及映射波数k
′
：
[0029][0030]
式(4)中：k为波数。
[0031]
在上述技术方案中，当空间偏导数选用空间二阶偏导数时，
[0032]
空间二阶偏导数的2n阶精度中心差分离散可以表示为：
[0033][0034]
式(5)中：p(x,t)为关于空间坐标x以及时间变量t的波场变量；pi表示t时刻位于空间离散位置iδx处的波场；δx为空间离散步长；cl为高阶差分算子系数；
[0035]
空间二阶偏导数的2n阶精度中心差分离散波数映射关系及映射波数k
′
：
[0036][0037]
式(6)中：k为波数。
[0038]
在上述技术方案中，在步骤三中，对震源波场进行空间频散正变换(fsdt计算)，包括如下步骤：
[0039]
步骤3.1：依据不同差分离散格式，代入相对应的映射波数k
′
，对参考震源波场p(x
′
,t0)进行离散傅里叶变换，得到频率域波场p(k,t0)，公式如下：
[0040][0041]
式(7)中：x为空间传播距离长度；
[0042]
步骤3.2：对频率域波场p(k,t0)进行快速傅里叶算法，反变换回时间域得到具有数值频散的波场p
′
(x,t0)，公式如下：
[0043][0044]
式(8)中：δx为空间离散步长；k为波数。
[0045]
在上述技术方案中，在步骤四中，数值频散波场与参考波场之间的归一化二范数误差设为频散误差，有如下表达式：
[0046][0047]
式(9)中：p(x
′
,t0)为参考波场；p
′
(x,t0)加完数值频散的波场；
[0048]
用数值频散波场与真实波场之间误差作为评判频散程度的依据，根据频散误差大小选取合适的误差阈值。
[0049]
在上述技术方案中，在步骤五中，筛选有限差分参数组合的具体方法为：根据步骤四中选取的误差阈值，筛选出误差小于或等于误差阈值时所对应的所有有限差分参数组合。
[0050]
在上述技术方案中，在步骤四中，阈值选取为0.01；
[0051]
在步骤五中，筛选有限差分参数组合的具体方法为：依据步骤四中选取的误差阈值0.01，筛选出误差小于或等于0.01时所对应的所有有限差分参数组合。
[0052]
在上述技术方案中，在步骤六中，设定目标准则，优化选取有限差分参数；
[0053]
目标准则包括计算量、存储量、flops、macs和flops等，依据目标准则的不同对有限差分参数组进行筛选，优化选取有限差分参数，具体方法为：
[0054]
步骤6.1：设定目标准则(例如设定计算量最小目标准则，其他目标准则类似)；
[0055]
步骤6.2：对模型进行网格划分，统计限差分参数组的目标准则相对应值(例如统计有限差分参数组相对应的计算量)；
[0056]
步骤6.3：对目标准则的相对应值进行统计，依据目标准则的相对应值优化选取出有限差分参数(例如对计算量进行统计，选取计算量最小时相对应的有限差分参数)。
[0057]
本发明由于采用以上技术方案，其具有的以下优点：
[0058]
(1)利用空间频散正变换(fsdt)方法可以高效率计算得到数值频散波场；
[0059]
(2)给定一个频散误差阈值，可以计算相对应的差分算子长度和波长采样点数等有限差分参数组合；
[0060]
(3)基于fsdt方法计算的数值频散误差，可以为选取最优的有限差分参数提供理论指导。
附图说明
[0061]
图1是本发明的整体流程示意图；
[0062]
图2为本发明实施例中计算量最小目标对应的有限差分参数组合示意图；
[0063]
图2中g表示最小波长采样点数。
具体实施方式
[0064]
下面结合附图详细说明本发明的实施情况，但它们并不构成对本发明的限定，仅作举例而已。同时通过说明使本发明的优点更加清楚和容易理解。
[0065]
本发明是优化选取出有限差分参数，是依据一定的目标准则而不是依据优化差分参数，本发明为有限差分参数优化选取提供理论指导。本发明给出了可以高效计算出波场空间方向数值频散的解析表达式，提出了空间频散正变换(forwardspacedispersiontransform，fsdt)方法理论，基于满足一定的目标准则(例如计算量最小、存储量最小、flops最小、macs最小、flops尽可能最大等)和空间频散误差准则条件下，可以优化的选取出有限差分参数，保证在一定数值模拟精度的条件下，使有限差分数值模拟的效率最高，节省计算成本，促进生产效益。
[0066]
如图1所示，本发明提供了有限差分数值模拟频散和参数优化选取方法，其包括以下步骤：
[0067]
1)根据真实地震记录提取出震源子波或者理论震源子波(例如ricker子波)；
[0068]
2)对空间偏导数进行有限差分离散，得到波数映射关系以及映射波数；
[0069]
具体如下：
[0070]
对空间偏导数进行有限差分离散，有限差分离散格式不同对应的波数映射关系也不同。这里仅给出了空间一阶、二阶偏导数的2n阶精度中心差分离散格式以及空间一阶偏导数的2n阶精度交错网格差分离散格式，其他差分格式具有类似推导。
[0071]
例如，空间一阶偏导数的2n阶精度中心差分离散可以表示为：
[0072][0073]
式(1)中：p(x,t)为关于空间坐标x以及时间变量t的波场变量；p
i l
表示t时刻位于空间离散位置(i l)δx处的波场；δx为空间离散步长；cl为高阶差分算子系数；
[0074]
例如，空间二阶偏导数的2n阶精度中心差分离散可以表示为：
[0075][0076]
例如，空间一阶偏导数的2n阶精度交错网格差分离散可以表示为：
[0077][0078]
式中，p(x,t)为关于空间坐标x以及时间变量t的波场变量；pi表示t时刻位于空间离散位置iδx处的波场；p
i l-1/2
表示t时刻位于空间离散位置(i l-1/2)δx处的波场；δx为空间离散步长；c
l
为高阶差分算子系数。对相关差分离散表达式进行傅里叶变换及简化，得到波数k的映射关系及映射波数k
′
：
[0079]
例如，空间一阶偏导数的2n阶精度中心差分离散，其波数k的映射关系及映射波数k
′
可以表示为：
[0080][0081]
式(2)中：k为波数；
[0082]
例如，空间二阶偏导数的2n阶精度中心差分离散，其正值波数k的映射关系及映射波数k
′
可以表示为：
[0083][0084]
式(6)中：k为波数；
[0085]
例如，空间一阶偏导数的2n阶精度交错网格差分离散，其波数k的映射关系及映射波数k
′
可以表示为：
[0086][0087]
式(4)中：k为波数；
[0088]
3)提出空间频散正变换(forward space dispersion transform，fsdt)方法理论，对震源波场进行fsdt计算；
[0089]
具体如下：
[0090]
设置没有数值频散的震源波场p(x
′
,t0)为参考波场，对参考波场进行空间频散正变换(forward space dispersion transform，简称fsdt)，可以表示为：
[0091][0092]
式中，k
′
为正映射波数，t0为时间变量，x为空间坐标变量，x为波场传播距离。fsdt方法具体实施过程为：
[0093]
3.1)代入映射波数k
′
，对参考波场p(x
′
,t0)进行离散傅里叶变换，得到频率域波场p(k,t0)，公式如下：
[0094][0095]
式(7)中：x为空间传播距离长度，
[0096]
3.2)对频率域波场p(k,t0)进行快速傅里叶算法，反变换回时间域得到具有数值频散的波场p
′
(x,t0)，公式如下：
[0097][0098]
式(8)中：δx为空间离散步长，k为波数；
[0099]
4)设定频散误差评定准则，选取合适误差阈值；
[0100]
具体如下：
[0101]
建议数值频散波场与参考波场之间的归一化二范数误差设为频散误差，并给出如下频散误差表达式：
[0102][0103]
式(9)中：p(x
′
,t0)为参考波场，p
′
(x,t0)加完数值频散的波场。
[0104]
建议用数值频散波场和真实波场之间误差作为评判频散程度的依据，根据频散误差大小选取合适的误差阈值，建议阈值选取为0.01。
[0105]
5)给定合适误差阈值，对不同组合的有限差分参数(例如：空间网格划分精度、有限差分算子长度、有限差分系数)进行筛选，得到给定误差阈值下的有限差分参数组合。
[0106]
6)设定目标准则(例如：计算量、存储量、flops、macs、flops等)，依据目标准则对不同的有限差分参数组进行筛选，优化选取有限差分参数，图2所示为计算量最小目标对应的有限差分参数组合，图中g表示最小波长采样点数。图2(a)为网格划分精度和有限差分算子长度的关系示意图；图2(b)为网格划分精度和计算量的关系示意图。
[0107]
具体如下：
[0108]
6.1)设定目标准则，例如设定计算量最小目标。
[0109]
6.2)对模型进行网格划分，统计有限差分参数组相对应的计算量。
[0110]
6.3)对计算量进行统计，依据计算量最小目标优化选取出有限差分参数。从图2中的(b)图可以看出，计算量曲线最小值点即为最优化的网格划分间距。从图2中的(a)图可以看出，最优化的网格划分间距对应的算子长度即为最优化有限差分算子阶数。
[0111]
综上所述，本发明提出的空间频散正变换方法具有解析表达式，可以高效率计算出波场的有限差分数值频散，依据频散误差阈值和目标准则，可以优化选取有限差分参数，进一步促进有限差分方法发展。
[0112]
实施例
[0113]
依据本发明提供的空间数值频散模拟方法和有限差分参数优化选取方法对某雷克子波的空间数值频散模拟和有限差分参数优化选取，选取震源子波为雷克子波时(其他震源参考波场类似)，有限差分算子阶数和网格划分精度关系曲线以及网格划分精度和计算量关系曲线如图2所示。图2中的(a)图表示有限差分算子阶数和网格划分精度关系曲线图；图2中的(b)图以及网格划分精度和计算量关系曲线图。从图2中的(b)图可以看出，计算量曲线最小值点即为最优化的网格划分间距。再从图2中的(a)图可以看出，最优化的网格划分间距对应的算子长度即为最优化有限差分算子阶数。依据不同的目标：例如存储量、flops、macs、flops等，均可以得到类似图2的关系曲线，从图2中可以很容易得到最优化的有限差分模拟参数(有限差分算子阶数为18阶，模型网格划分精度g约为2.6，其图2中g表示雷克子波2.5倍主频波长对应的采样点数)。
[0114]
上述各实施例仅用于说明本发明，各个步骤都是可以有所变化的，在本发明技术方案的基础上，凡根据本发明原理对个别步骤进行的改进和等同变换，或者依据本发明原理对个别步骤进行的改进和等同变换的算法、软件以及应用，均不应排除在本发明的保护范围之外。
[0115]
其它未说明的部分均属于现有技术。