一种任意形态不连续面自适应表达的方法与流程

1.本发明属于地震勘探技术领域，尤其涉及一种任意形态不连续面自适应表达的方法。


背景技术：

2.在地震学中，通常把地震传播过程中突然发生变化的地方成为不连续面。根据不连续面两侧介质物理的物理性质(固体，液体和气体)，可以将不连续面划分为以下四类：1)固体/固体不连续面；2)固体/液体不连续面；3)固体自由表面(陆地表面)，也称固体/真空不连续面；4)液体自由表面(海平面)，也称液体/真空不连续面。这类不连续面附近均会发生复杂的波转换现象和强能量的界面波生成。例如，固体/液体不连续面往往用于描述海底界面或者地球内外核边界，在其附近会产生一种沿该界面传播的scholte wave，以及发生p
‑
s波转换，会产生强振幅面波(rayleigh wave或love wave)，其频散特征在近地表面地震勘探中有重要的意义和用途；液体自由表面用于描述海洋表面，是海洋勘探中多次波产生的原因，其可以被视为噪音进行除去(多次波去噪技术)，也可以用于照明角度和区域补偿(多次波成像技术)。
3.在复杂的介质地震波数值模拟中，往往会遇到以上不连续面的组合，对于这些不连续面的正确数值建模表达会直接影响数值模拟的准确性和精度，特别是对于地表观测地震数据和海底地震仪观测地震数据，其检波器和地震仪直接位于不连续面上或附附近，因此受不连续面的影响更大。能否正确的对不连续面附近的波长进行数值模拟会直接对这类地震数据中提取和分析信息产生影响。
4.此外，这类不连续面还会破坏一些数值离散方法的完整性，例如有限差分法，伪谱法等基于偏微分方程强解形式的算法。目前工业界和科研界最常用的地震数值模拟技术是有限差分法，其适用范围为参数均匀或光滑变化介质，直接应用于含强阻抗对比的不连续面时会导致计算不稳定，数值假象产生和计算精度降低，进而影响最终模拟结果的准确性。
5.目前，一般采用直接离散化边界条件，以差分形式来近似边界条件中的控件导数，以固体自由表面的实现为例，有中心差分显示近似格式(alterman和karal，1968)、单边显示近似格式(vidale和clayton，1986)以及复合近似等格式(ilan等，1975；ilan和loewenthal，1976；lan和zhang，2011)；另外，在一阶速度
‑
应力方程表达形式下对应力张量和速度矢量分别操作，以固体自由表面的实现为例，有真空格式(zahradnik等，1993；bohlen和saenger，2006)、应力镜像法(levander，1988；graves，1996；kristek等，2002)和一些介质平均方法(mittet，2002；xu等，2007)；更有，基于mimetic算子的离散方法，例如dela puente等(2014)，shragge和tapley(2017)，konuk和shragge(2020)，shthi等(2021)；除以上三种方法，还有有限元类方法，包含谱yuan法(komatitsch和tromp，1999；komatitsch等，1999；tromp等，2008；peter等，2011)，这类基于偏微分弱解形式的数值模拟方法，可参考现有的普元发开源程序包：specfem2d和specfem3d。
6.但是，现有技术存在以下问题：1)应用于高泊松比介质的正演模拟会产生计算不
稳定的现象，并且最高只能到达二阶的计算精度。在截面起伏情况下其计算稳定性还会受到截面最大曲率的影响，曲率越大越容易不稳定；2)计算精度低，且受泊松比变化影响；在截面起伏情况下，截面两侧物理场所需的关于界面的对称操作变得及其复杂；3)算法实现复杂，需额外开辟内存控件来存储集合拓扑变量，且计算稳定性还会收到截面最大曲率的影响，曲率越大越容易不稳定；4)相比较于基于有限差分法的数值模拟技术，算法理论和实现复杂，计算量大，内存消耗费多。


技术实现要素：

7.针对上述背景技术中提到的技术缺陷，为了解决上述题述问题，本发明提供了一种任意形态不连续面自适应表达的方法，其特征在于，包括如下步骤：
8.s1：导入初始正演模型，包括纵横波的速度和密度，设置震源子波的类型和主频的参数，并根据各模型的参数设置空间步长和时间步长；
9.s2：根据模型参数选取时间步长和空间步长，设置差分的阶数和震源子波；
10.s3：当模型表面起伏时，对模型的自由表面进行阶梯化离散，将非规则的起伏地表变成规则的阶梯状网格；
11.s4：运用介质参数修正法对自由表面上非水平的点进行本构关系和密度的修正；
12.s5：利用修正后的本构关系来替换自由表面处原来的本构关系，在不对原始的有限差分法代码进行修改的情况下，即可引入自由表面的影响；
13.s6：对于各向异性介质，将修正后的本构关系带入到有限差分方法的代码中，即可在各向异性介质的情况下引入自由表面的影响。
14.进一步的，所述自由表面包括液体自由表面和固体自由表面。
15.进一步的，所述不连续面包括固体/液体不连续面和液体/真空不连续面。
16.进一步的，对于各异性介质，需要导入各向异性参数。
17.进一步的，所述异性介质适用于固体自由表面部分。
18.本发明所述通过介质参数修正实现任意形态不连续面自适应表达的方法，可直接继承与应用于现有的工业界和学术界主流的速度
‑
应力交错网格有限分数值程序，精确贴合不连续面的空间位置，没有传统方法中的版网格偏移误差；通过简单的修正不连续面附近网格点的介质参数来实现，且该修正过程只需一次，计算高效、操作简单；可应用高阶有限差分算子，实现与基于偏微分方程弱解类精确实现方法相当的计算精度；适用于各向异性介质附近的不连续面，且具有泊松比自适应的特点。
附图说明
19.为了更清楚地说明本发明实施例中的技术方案，下面对实施例描述中所需要使用的附图作一个简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
20.图1为本发明实施例的介质参数修正实现任意形态不连续面自适应表达的方法的流程图；
21.图2为本发明实施例的自由表面离散化示意图；
22.图3为本发明实施例的固体与空气的位置关系示意图，图a)指空气在自由表面左侧，图b)指空气在自由表面前面，图c)空气在自由表面右侧，图d)指空气在自由表面后面，图e)指空气在自由表面上面。
具体实施方式
23.为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
24.本发明提供了一种任意形态不连续面自适应表达的方法，其特征在于，包括如下步骤：
25.s1：导入初始正演模型，包括纵横波的速度和密度，设置震源子波的类型和主频的参数，并根据各模型的参数设置空间步长和时间步长；
26.s2：根据模型参数选取时间步长和空间步长，设置差分的阶数和震源子波；
27.s3：当模型表面起伏时，对模型的自由表面进行阶梯化离散，将非规则的起伏地表变成规则的阶梯状网格；
28.s4：运用介质参数修正法对自由表面上非水平的点进行本构关系和密度的修正；
29.s5：利用修正后的本构关系来替换自由表面处原来的本构关系，在不对原始的有限差分法代码进行修改的情况下，即可引入自由表面的影响，所述影响包括但不限于自由表面引起的地震波的反射；
30.s6：对于各向异性介质，将修正后的本构关系带入到有限差分方法的代码中，即可在各向异性介质的情况下引入自由表面的影响。
31.进一步的，所述自由表面包括液体自由表面和固体自由表面；所述不连续面包括固体/液体不连续面和液体/真空不连续面；所述异性介质适用于固体自由表面部分。
32.具体地，对固体/液体不连续面假设固体介质1对应的本构关系为固体介质2对应的本构关系为根据平均介质思想，可将两个连接的介质等效成一个等效平均介质，其本构关系为同理，对于固体/液体不连续面，假设液体介质1对应的本构关系为液体介质2对应的本构关系为根据平均介质思想，可将两个连接的介质等效成一个等效平均介质，其本构关系为
33.上述式中：τ指应力张量矩阵；m1指弹性介质1，m2指弹性介质2，a代表平均介质；e指等效的弹性系数矩阵；λ1、λ2和μ1、μ2是拉梅常数；ε为应变张量矩阵。
34.进一步的，对于固体自由表面，例如陆地表面，即固体/真空不连续面，假设上层真空介质的本构关系为其中根据λ1→
0，μ1＝0，下层固体介质的本构关系为则等效平均介质的本构关系为λ1→
0，μ1＝0；
35.进一步的，对于地下各向异性介质情况，即固体/真空不连续面，假设上层真空介质的本构关系为其中对于下层各向异性介质，其本构关系为则等效平均介质的本构关系为其中
36.上述公式中c
11，
c
13，
c
33
，c
44
为各向异性介质的弹性参数。
37.上述式中：τ指应力张量矩阵；c代表弹性系数；m1指弹性介质1，m2指弹性介质2，a代表平均介质；e指等效的弹性系数矩阵；λ1、λ2和μ1、μ2是拉梅常数；ε为应变张量矩阵。
38.具体的，在步骤s1中，对于各向异性介质，需要导入各向异性参数c
11，
c
13，
c
33
，c
44
。同时需要确定差分的阶数、震源类型和震源主频等参数，根据模型参数设置空间步长和时间步长。
39.进一步的，当模型的地表面起伏时，需要对模型的自由面进行阶梯化离散，离散如图1所述，将非规则的欺负地面变成规则的阶梯状网格，如图2，图2中0l、0r指网格点的外角点，il、ir代表网格点的内角点，图2显示了起伏地表的实现方式；对自由表面网格点进行划分，并对不同网格点根据表1进行本构关系修正。
40.在步骤s4中，运用介质参数法对自由表面上非水平的点进行本构关系修正，修正关系如下表所示，其中vl代表自由表面左边为空气的情况，vr代表自由表面右边为空气的情况，vb代表自由表面后边为空气的情况。h代表自由表面上方为空气的情况，具体可参考图3，图3中a指空气，s指地下介质。
41.表1：起伏地表自由表面处本构关系
[0042][0043]
上述表中：τ
ij
指应力张量的第ij分量；ε
ij
指应变张量的第ij分量；λ和μ指拉梅常数；ρ指介质的密度；α＝(2μ(λ+μ)/(λ+2μ))，β＝(μλ)/(λ+2μ)。
[0044]
进一步的，利用修正后的本构关系，只需更改模型的第一层参数，即自由表面处的参数，这样在不对原始的有限差分法代码进行修改的情况下就可以引入自由表面的影响。在自由表面处的本构关系如下所示：
[0045][0046]
[0047][0048][0049][0050]
τ
zz|z
＝0＝0，
[0051]
τ
xy|z＝0
＝τ
yx|z＝0
＝με
yx
＝με
xy
，
[0052]
τ
yz|z＝0
＝τ
zy|z＝0
＝0，
[0053]
τ
xz|z＝
0＝τ
zx|z＝0
＝0，
[0054]
u
i
是位移的第i个分量，τ
ij
指应力张量的第ij分量；ε
ij
指应变张量的第ij分量；λ和μ指拉梅常数；ρ指介质的密度。
[0055]
进一步的，对于各向异性介质，以vti介质为例，vti为一种各向异性介质，其修正后的本构关系如下式所示：
[0056][0057]
修正后的本构关系带入到位移应力方程中可得：
[0058]
[0059][0060]
由此即可在各向异性介质的情况下引入自由表面的影响。
[0061]
本发明所述通过介质参数修正实现任意形态不连续面自适应表达的方法，可直接继承与应用于现有的工业界和学术界主流的速度
‑
应力交错网格有限分数值程序，精确贴合不连续面的空间位置，没有传统方法中的版网格偏移误差；通过简单的修正不连续面附近网格点的介质参数来实现，且该修正过程只需一次，计算高效、操作简单；可应用高阶有限差分算子，实现与基于偏微分方程弱解类精确实现方法相当的计算精度；适用于各向异性介质附近的不连续面，且具有泊松比自适应的特点。
[0062]
尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例和落入本发明范围的所有变更和修改。
[0063]
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。