一种基于分数阶拉氏算子的各向异性衰减介质模拟方法

1.本技术涉及油气探测领域，具体涉及一种基于分数阶拉氏算子的各向异性衰减介质模拟方法。


背景技术：

2.在油气探测开发中，地震勘探技术发挥着极其重要的作用。目前，发展针对复杂介质的地震波场数值模拟已经成为当今工业界的共识，近年来发展的分数阶解耦(dfl)粘滞波动方程具有诸多独有优势。然而，dfl粘滞波动方程的数值求解往往面临着难题。该方程包含的阶数随空间变化的拉普拉斯算子给非均匀介质的波场模拟带来了困难。目前，虽然已经存在dfl形式的各向异性粘滞弹性波和声波方程，但其中包含的混合域拉普拉斯算子需要额外的处理才能适用于非均匀介质的模拟。
3.专利文献(cn113341455a)公开了一种粘滞各向介质地震波数值模拟方法，该技术方案虽然能够同时准确描述速度和衰减各向异性、同时又便于数值计算方程，但是其只考虑了纵波的衰减和速度的各向异性，存在不符合实际地下情况的问题，会使得地震模拟不准确，因为地震波在地下同时以纵波和横波传递。


技术实现要素：

4.针对上述技术问题，本技术提供一种基于分数阶拉氏算子的各向异性衰减介质模拟方法，同时考虑了纵波和横波的各向异性，能够更好的贴合地下实际情况，使得地震模拟更准确。
5.本发明实施例采用的技术方案为：
6.本发明实施例提供一种基于分数阶拉氏算子的各向异性衰减介质模拟方法，所述方法包括如下步骤：
7.s100，基于预设的频散关系构建弹性波的频散关系，所述弹性波包括p 波和s波；
8.s110，基于构建的弹性波的频散关系得到弹性衰减介质的复模量；
9.s120，获取vti衰减介质的本构关系在频率-波数域的表达式；
10.s130，将得到的弹性衰减介质的复模量扩展至vti衰减各向异性，将 vti衰减介质的本构关系在频率-波数域的表达式转换为在时间-空间域的表达式；
11.s140，基于vti衰减介质的本构关系在时间-空间域的表达式、几何方程和运动平衡方程，得到vti介质dfl粘滞弹性方程；
12.确定地下介质参数，并将所述地下介质参数和地震波参数代入到所述 vti介质dfl粘滞弹性方程中，以计算获得地震波波场模拟数值；
13.其中，所述vti介质dfl粘滞弹性方程为：
[0014][0015]
其中，t表示时间，v
x
和vz分别为地下介质质点沿x方向和z方向的振动速度分量，ρ为地下介质的密度，σ
xx
、σ
zz
和σ
xz
分别为应力张量的水平分量、纵向分量和剪切分量；
[0016]
其中，
[0017]q11
和q
33
分别为p波在横向和垂直方向上的品质因子， q
55
为s波的品质因子，
[0018]vp
为p波沿对称轴方向的速度；vs是s波的速度，ε和δ为速度各向异性参数，δq为衰减各向异性参数；为p波沿水平对阵轴的速度，v2＝v
p
，为p波沿垂直对阵轴的速度；ω0为参考角频率。
[0019]
本发明至少具有以下技术效果：本发明提供的基于分数阶拉氏算子的各向异性衰减介质模拟方法，同时考虑了纵波和横波的各向异性，能够更好的贴合地下实际情况，使得能更准确地刻画地震波在地下传播过程中的衰减机制。
附图说明
[0020]
为了更清楚地说明本技术实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本技术的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0021]
图1为本技术实施例提供的一种基于分数阶拉氏算子的各向异性衰减介质模拟方法的流程图；
[0022]
图2(a)和图2(b)为通过与解析解对比来验证本技术实施例提供的 vti介质dfl粘滞弹性方程的模拟精度的示意图；
[0023]
图3为验证本技术实施例提供的vti介质dfl粘滞弹性方程的解耦特性的示意图；
[0024]
图4a～c为检验本技术实施例提供的vti介质dfl粘滞弹性方程处理q 值和速度剧烈变化介质的能力的三个层状模型；
[0025]
图5为本技术实施例提供的谱比法估算q值的观测系统图；
[0026]
图6为本技术实施例提供的谱比法所用观测系统记录的共炮点道集示意图；
[0027]
图7为本技术实施例提供的品质因子的倒数对比图；
[0028]
图8a至图8d为本技术实施例提供的2007bp tti模型示意图；
[0029]
图9a为本技术实施例提供的各向异性弹性介质中的波场快照；
[0030]
图9b和图9c分别为各向异性dfl粘滞弹性波方程和各向异性时间分数阶粘滞弹性波方程的模拟结果示意图；
[0031]
图9d为图9b和图9c的差值；
[0032]
图10a为本技术实施例提供的各向异性弹性介质中的地震记录示意图；
[0033]
图10b和图10c分别为各向异性dfl粘滞弹性波方程和各向异性时间分数阶粘滞弹性波方程的模拟结果示意图；
[0034]
图10d为图10b和图10c的差值；
[0035]
图11a为索尔顿槽模型的位置图；
[0036]
图11b为重新采样的p波速度模型示意图；
[0037]
图11c和图11d分别为索尔顿槽模型不同纬度和不同经度的剖面示意图；
[0038]
图12a至图12f为本技术实施例提供的dfl粘弹性波动方程计算不同时刻的波场快照图；
[0039]
图13a～c分别为采用本技术实施例提供的粘弹性波动方程、现有的时间分数阶粘弹性波动方程弹性波动方程计算三维共炮点道集示意图；
[0040]
图14为本技术实施例提供的粘弹性波动方程与v
x
分量、vy分量和vz分量的地震记录中视频包络目标函数、时间分数阶粘弹性波动方程和时频相位误差进行跟踪比较的示意图。
具体实施方式
[0041]
下面将结合本技术实施例中的附图，对本技术实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本技术一部分实施例，而不是全部的实施例。基于本技术中的实施例，本领域技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本技术保护的范围。
[0042]
如图1所示，本发明实施例提供的基于分数阶拉氏算子的各向异性衰减介质模拟方法包括如下步骤：
[0043]
s100，基于预设的频散关系构建弹性波的频散关系，所述弹性波包括p 波和s波。
[0044]
在本发明实施例中，所述预设的频散关系为常q模型的近似频散关系。常q模型的近似频散关系表示频率和速度的变化关系，其公式为：
[0045]
ω2＝d1k+d2k2+d3k3+d4(iω)k+d5(iω)k2ꢀꢀꢀ
(1)
[0046]
(d1,d2,d3,d4,d5)＝(-γvω0,v2,γv3ω
0-1
,πγv,πγ2v2ω
0-1
)
ꢀꢀꢀ
(2)
[0047]
其中，k为波数，ω为角频率，ω0为参考角频率，γ为空变分数阶阶数，γ＝arctanq-1
/π，q为弹性波的品质因子；v为参考速度，v＝v0cos(πγ/2)，v0为在参考频率ω0处的参考相速度。
[0048]
弹性波同时包含p波和s波分量，而二者都满足上述近似的频散关系，因此，构建的弹性波的常q模型的频散关系为：
[0049][0050]
其中，v
α
为弹性波的速度，γ
α
＝arctanq
α-1
/π，q
α
为弹性波的品质因子，α＝p或s。
[0051]
s110，基于构建的弹性波的频散关系得到弹性衰减介质的复模量。
[0052]
在本发明实施例中，所述弹性衰减介质的复模量可为：
[0053][0054]
s120，获取vti衰减介质的本构关系在频率-波数域的表达式。
[0055]
在本技术实施例中，考虑较为简单的各向异性介质：具有垂直对称轴的横向各向同性介质(vti)。其应力应变关系是由与介质弹性性质密切相关的弹性刚度矩阵决定的，对于vti弹性介质，其vogit形式的弹性刚度矩阵的表达式为
[0056][0057]
弹性刚度系数为：
[0058][0059]
其中ρ表示地下介质的密度；v
p
表示p波沿对称轴方向的速度；vs是s波的速度，ε和δ为速度各向异性参数。
[0060]
对于vti衰减介质，其q矩阵可以表示为
[0061][0062]
其中q
11
和q
33
分别对应p波在横向和垂向上的品质因子，而q
55
对应s 波的品质因子，q
12
＝q
13
,q
55
＝q
66
。此外，使用类似于thomsen系数的衰减各向异性参数εq和δq，以此来建立起q
11
，q
33
和q
13
之间的关系，具体地，
[0063][0064][0065]
与各向异性弹性介质的刚度系数为实数不同，各向异性衰减介质对应的是复刚度系数(或复模量，记为)，它可以通过刚度矩阵由如下关系转换得到
[0066][0067]
vti衰减介质的本构关系在频率-波数域的表达式为
[0068][0069]
[0070][0071]
其中，和分别为应力张量的水平分量、纵向分量和剪切分量在频率中的表达式，和分别为应变张量的水平分量、纵向分量和剪切分量在频率中的表达式，和分别为各向异性衰减介质的复模量。
[0072]
s130，将得到的弹性衰减介质的复模量扩展至vti衰减各向异性，将 vti衰减介质的本构关系在频率-波数域的表达式转换为在时间-空间域的表达式。
[0073]
该步骤具体包括：
[0074]
s131，将式(3)分解为p波和s波分量，得到对应的和其中，
[0075][0076][0077]
s132，因为式(14)只与v1有关，因此，根据式(4)可推导出式(14) 中的的近似表达式，具体地，将v1替换式(4)中的v
α
，得到式(14)中的的近似表达式为：
[0078][0079]
s133，因为式(15)只与vs有关，因此，根据式(4)可推导出式(15) 中的和的近似表达式，具体地，将vs替换式(1)中的v
α
，得到式(15) 中的和的近似表达式为：
[0080][0081][0082]
s134，将式(16)代入式(14)以及将式(17)和式(18)代入式(15)，得到式(11)在时间-空间域的表达式：
[0083][0084]
s135，对式(12)和(13)分别执行与s131～s134类似的操作，即分别分解为p波和s波分量，然后基于式(4)分别得到式(12)和(13)在时间
ꢀ‑
空间域的表达式：
[0085][0086]
σ
xz
＝l
55exz
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(21)
[0087]
其中，e
xx
、e
zz
和e
xz
分别为应变张量的水平分量、纵向分量和剪切分量。
[0088]
s140，基于vti衰减介质的本构关系在时间-空间域的表达式、几何方程和运动平衡方程，得到vti介质dfl粘滞弹性方程。
[0089]
在本发明实施例中，所述几何方程和运动平衡方程为：
[0090][0091]
其中，i和j为x或者z。
[0092]
在本发明实施例中，结合式(19)至(22)得到的基于新频散关系的 vti介质dfl粘滞弹性方程为：
[0093][0094]
其中，t表示时间，v
x
和vz分别为地下介质质点沿x方向和z方向的振动速度分量，σ
xx
、σ
zz
和σ
xz
分别为应力张量的水平分量、纵向分量和剪切分量；
[0095]
其中，
[0096]vp
为p波沿对称轴方向的速度；vs是s波的速度，ε和δ为速度各
向异性参数，为p波沿水平对阵轴的速度，v2＝v
p
，为p波沿垂直对阵轴的速度。
[0097]
s150，确定地下介质参数，并将所述地下介质参数代入到所述vti介质 dfl粘滞弹性方程中，以计算获得地震波波场模拟数值。
[0098]
在得到vti介质dfl粘滞弹性方程的基础上，通过确定地下介质参数包括地下介质密度、p波和s波品质因子、p波速度和s波速度、衰减各向异性参数、速度各向异性参数等，并将上述参数代入到vti介质dfl粘滞弹性方程，从而实现了地震波波场模拟，得到地震波波场模拟数值。
[0099]
(测试例1)
[0100]
在该测试例中，首先通过与解析解进行对比验证上述式(23)示出的 dfl粘滞弹性波方程的精度，其中，在均匀介质中的解析解可以由格林函数推导得到。所使用的均匀介质模型大小的为512
×
512个网格点，空间采样间隔为10m。设置在均匀介质中心的雷克子波震源主频为25hz，参考频率与震源一致的参考速度分别为c
0p
＝3000m/s和c
0s
＝2000m/s，介质密度为2200 kg/m3。记录(3160m，3160m)处的v
x
分量，图2a部分和图2b部分分别对应完全弹性介质和衰减介质(q
p
＝32，qs＝20)。在图中，粗线、虚线以及细黑线分别对应解析解、数值解以及二者残差。可以发现，无论对于衰减介质还是完全弹性介质，所提出新的方程都能非常好地吻合解析解，证明了dfl粘滞弹性波方程的准确性。
[0101]
(测试例2)
[0102]
在该测试例中，进一步验证dfl粘滞弹性波方程的解耦特性。在图3中，四个波场快照分别对应于完全弹性波方程、衰减主控方程、频散主控方程和所提出的dfl粘滞弹性波方程。可以观察到，与弹性波相比，衰减主控方程很好地描述了振幅能量的衰减，频散主控方程刻画了地层衰减造成的速度频散(可观察到相位的滞后)，而dfl粘滞弹性波方程模拟的波场则同时具有这两种效应。
[0103]
(测试例3)
[0104]
在该测试例中，进一步验证dfl粘滞弹性波方程处理非均匀衰减介质的能力。在该测试例中，设计三个层状模型来检验其处理q值和速度剧烈变化介质的能力。如图4(a)所示，模型1中的速度为层状介质而q值均匀，模型 2中的速度均匀而q值为层状，模型3中的速度和q均为层状变化。该组模型大小为800
×
800个网格点，空间采样间隔为10m，水平界面位于深度 4400m处，采用的时间步长为1ms，设置在介质中心的雷克子波震源主频为 25hz。然而，在q对比明显的情况下(模型2和模型3)，方案之间的快照残差比参考快照更显著，这是相位畸变叠加效应的结果(见图4c)。图4b中的第三列只存在细微的残差，图4c中的两条区线匹配良好，这些都证明了本发明提供的dfl粘滞弹性波方程处理非均匀衰减介质的良好能力。
[0105]
(测试例4)
[0106]
在该测试例中，为验证所提出的vti-dfl粘滞弹性方程描述衰减随方向变化的准确性，进行了一组各向异性均质模型试验。首先求解提出的vti
‑ꢀ
dfl粘滞弹性方程，得到相应的快照。一方面，利用谱比法可以估计各方向的q值；另一方面，通过求解平面波方程的christoffel方程，可以得到理论 q值。通过比较理论q值与估计q值的相容性，评价本发明提出的vti-dfl 粘滞弹性方程的准确性。图5为仿真几何图形，其中星形为主频为25hz的雷
克子波，两个由三角形组成的同心圆为接收器，每个角度间隔600m分布。测试所使用的是一组均匀各向异性介质模型，其衰减各向异性参数如表1所示。模型大小为400
×
400个网格点，空间采样间隔10m，时间步长1ms。其余参数为v
p
＝3.0km/s，vs＝1.5km/s，thomsen速度各向异性参数为ε＝0.16，δ＝0.08，介质密度2000kg/m3。为了定量地表述介质衰减随方向变化的强烈程度，定义烈度系数strength(q
p,s
)＝1-min(q
p,s
)/max(q
p,s
)。
[0107]
表1一组衰减各向异性参数变化的弹性衰减介质模型
[0108][0109][0110]
利用谱比法所用观测系统记录的共炮点道集可入6所示，其中，在图6 中，左侧为内侧同心圆检波点，右侧为外侧同心圆检波点。得到的品质因子的倒数对比可如图7所示，其中，圆圈表示谱比法估算的q值，黑色实线为 christoffel方程得到的理论q值。
[0111]
(测试例5)
[0112]
在该测试例中，以一个复杂模型验证所提出的vti-dfl粘滞弹性方程的模拟精度，所采用模型为2007bp tti模型(不考虑倾斜角度)。图8a至图8d 所示分别为p波速度、p波品质因子q
p
，thomsen速度各向异性参数δ和ε。 s波速度和品质因子qs分别通过vs＝v
p
/1.73，qs＝q
p
/1.2转换得到，衰减各向异性参数为δq＝2δ，εq＝3ε。该模型大小为500
×
270个网格点，空间采样间隔10m，时间采样间隔1ms，模拟传播时间为2.5s，所采用的激发震源是主频为25hz的雷克子波，震源位置为(2500m,30m)，检波点位置位于沿水平方向50m深度的每个网格位置。为了评估模拟精度，我们从现有的vti-tf 粘弹性波动方程(zhu,2017)中计算了一个参考解，该方程能够描述与方向相关的衰减，但需要巨大的计算和内存开销。图9a所示为各向异性弹性介质中的波场快照，图9b由vti-dfl粘滞弹性方程模拟得到的波场快照，图9c由 vti介质时间分数阶粘滞弹性波方程模拟得到的波场快照，该方程在向各向异性拓展时不需要近似，所以可以很好地描述地层衰减随方向的变化，但是由于其中包含的时间分数阶算子在数值求解过程中存在计算和存储的问题，因此不适用于大规模的数值模拟；与该方程模拟结果对比，即可验证所提出 vti介质dfl粘滞方程的模拟精度。图9d为图9b和图9c相减的结果。图 10a至图10d分别为与图9a至图9d对应的共炮点道集地震记录。分析对比后可以发现，与非衰减介质相比，地层衰减效应会明显使振幅能量减弱；此外图9d和10d中微弱的差值表明，所提出的vti-dfl粘滞弹性方程对于复杂模型依然具有很高的数值模拟精度。
[0113]
(测试例6)
[0114]
在该测试例中，为了证明所提出的dfl粘弹性波动方程的通用性和可行性，在南加州索尔顿海槽的三维真实模型中模拟了地震波的传播。索尔顿海槽的定位图如图11a所示，最初的p波速度模型是由(ajala，2019)提供的，其中包含深度从0到9公里内的19个网格，经度115.7
°
w到116.7
°
w的 101网格，纬度从33.3
°
n 34.2
°
n内的91个网格。由于原始数据不是均匀的，进行线性插值和重采样可以得到一个新的p波速度模型，包含505
×
455
ꢀ×
48网格点，网格间距均匀为200m(见图11b)。图11c和11d分别显示了不同纬度和经度的一些剖面。s波速度由vs＝v
p
/1.73(li,1993)的速度模型得到， q模型由q
p
＝qs＝3.516(v
p
/1000)
22
(li,1993)的函数建立。在(50000m,46000m, 1600m)处放置一个主频为5hz的雷克子波震源，接收器位于地表x-y平面的每一个网格点。时间为15s，时间步长为5.0ms。图12a至图12f显示了采用所提出的dfl粘弹性波动方程模拟的分量在不同时间的快照。特别要注意的是，可以观察到一些能量被困在椭圆标记的区域，地震波在这个区域的传播速度也比周围地区慢。这种现象出现的原因可以从索尔顿海槽模型中找到，从图11b可以看出，该区浅部存在一些低速、高衰减的沉积物。为了评价所提出的索尔顿海槽模型dfl粘滞弹性方程的准确性，采用时间分数阶粘弹性波动方程作为参考解。由于计算时间分数阶粘弹性波动方程需要大量的计算机内存，为了使模拟可行，将模型重新采样到252
×
227
×
24网格点，网格间距均为20m。震源主频提高到25hz，时间步长减小到1.0ms。图13显示了 3d共炮点道集，其中水平切片显示了1.1s处的快照，两个垂直剖面显示了x 和y方向上的两个共炮点道集。图13a-c分别对应本发明提出的dfl粘弹性波动方程、时间分数阶粘弹性波动方程和弹性波动方程。图13a和13b中的共炮点道集看起来非常相似，由于能量衰减，它们的振幅明显弱于弹性情况 (图13c)。为了更好的对比，图14中进一步显示了(1500m,1500m,20m)的地震记录。从左到右的列对应的是v
x
，vy和vz的分量。上下两行分别表示时频包络目标函数(tfem)和时频相位误差(tfpm,kristekova et al.，2006，2009)。三个分量(图14中间行)的地震记录匹配良好，tfem和tfpm在2％以内，这表明所提出的dfl粘弹性波动方程生成的地震记录与参考解几乎相同。
[0115]
综上，本技术实施例提供的方法提供了基于分数阶拉普拉斯算子的声波方程，可以同时准确刻画地震波速度和衰减随方位变化的特性，而且不包含随空间变化的分数阶拉普拉斯算子，因此可天然地准确计算地震波在非均匀介质中的传播，通过对声波方程的简化处理，可以实现各向异性介质声波方程的降维处理(比如三维到二维)，从而可以广泛应用于二维(2d)和三维(3d)复杂衰减各向异性介质中声波波场模拟、吸收衰减补偿逆时偏移成像及全波形反演。并且，由于同时考虑了纵波和横波的各向异性，能够更好的贴合地下实际情况，使得能更准确地刻画地震波在地下传播过程中的衰减机制。
[0116]
本技术的实施例还提供了一种非瞬时性计算机可读存储介质，该存储介质可设置于电子设备之中以保存用于实现方法实施例中一种方法相关的至少一条指令或至少一段程序，该至少一条指令或该至少一段程序由该处理器加载并执行以实现上述实施例提供的方法。
[0117]
本技术的实施例还提供了一种电子设备，包括处理器和前述的非瞬时性计算机可读存储介质。
[0118]
本技术的实施例还提供一种计算机程序产品，其包括程序代码，当所述程序产品
在电子设备上运行时，所述程序代码用于使该电子设备执行本说明书上述描述的根据本技术各种示例性实施方式的方法中的步骤。
[0119]
虽然已经通过示例对本技术的一些特定实施例进行了详细说明，但是本领域的技术人员应该理解，以上示例仅是为了进行说明，而不是为了限制本技术的范围。本领域的技术人员还应理解，可以对实施例进行多种修改而不脱离本技术的范围和精神。本技术公开的范围由所附权利要求来限定。