一种基于节点间断有限元方法的波场正演模拟方法及系统

本发明属于勘探地球物理研究，具体涉及一种基于节点间断有限元方法的波场正演模拟方法及系统。

背景技术：

1、复杂地下模型地震波传播的数值模拟是全波场反演的重要基础，在地球物理中的地位日益重要。同时，平滑变化的物质在地质勘探中是非常常见的，例如沉积环境。然而，由于物性参数的多变性，在这种介质中很难得到高精度的数值解。

2、在过去的几十年里，在合成地震记录的计算中采用了各种各样的数值方法，如有限差分法(fdms)、有限单元法(fems)、谱元法(sems)。有限差分法由于其简单的数值格式和较高的计算效率，是地球物理中应用最广泛的数值方法。然而，fdms的几何灵活性较差，在包含特定界面包括声-弹耦合界面、弯曲界面和自由表面的介质中所得到解的精确度有限。有限元方法使用非结构性网格对计算域进行划分，该类网格能够清晰准确地表达介质中存在的各类界面，具有较高的几何灵活性，并且用低阶多项式逼近单元中的解。fems的主要缺点是，在每个时间步均需要对庞大的全局质量矩阵求逆，这极大地降低了方法的计算效率。谱元法是有限单元法的拓展，该方法在单元内使用基于gauss-lobatto-legendre(gll)插值点的高阶多项式对解进行近似。gll插值点的使用使得全局质量矩阵变为严格的对角矩阵，从而降低了质量矩阵求逆的计算难度，提高了计算效率。然而，由于采用了gll点，使得该方法更倾向于采用四边形或六面体单元进行网格划分，该类网格在模拟复杂地质构造时存在一定的生成困难。

3、间断有限元(dg)方法由于其局部性地求解过程并允许单元边界处解的不连续性，可以作为上述方法的一个可行的替代方案。dg方法最初由reed和hill提出，由hesthaven和warburton展开了进一步的数学研究。与有限体积法类似，该方法引入了数值流通量处理单元边界处的不连续性。kaser和dumbser最先将dg方法应用于地震波场正演，他们基于非结构网格提出了一种具有黎曼通量的任意高阶dg方法(ader-dg)。该方法后续被应用于粘弹性和各向异性介质。zhang等人通过使用godunov数值通量显式地处理单元界面，将ader-dg方法推广到多物理耦合介质中。常规的间断有限元法的缺点是物性参数在每一个单元内部必须都是恒定的。当单元尺寸接近波长时，连续变化的物性参数会导致模拟结果产生误差。节点间断有限元(ndg)方法在dg方法的基础上，在网格单元内部引入插值锚点，并使用基于锚点的拉格朗日插值多项式对局部解进行近似。在复杂模型中，单元内部物性参数不是恒定的情况下，节点间断有限元法可以得到更准确的结果。

技术实现思路

1、本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种基于节点间断有限元方法的波场正演模拟方法及系统，基于单元内部锚点处的物性参数对参数矩阵进行近似，实现对连续变化介质中地震波传播的计算和研究，并分析方法的精确性与可行性，用于解决间断有限元方法在处理连续变化介质模型时精度低、易出错的技术问题。

2、本发明采用以下技术方案：

3、一种基于节点间断有限元方法的波场正演模拟方法，包括以下步骤：

4、s1、对目标连续变化介质模型进行非结构网格划分，导出网格信息并关联相关物性参数，得到目标连续变化介质模型参数；

5、s2、确定波动方程；

6、s3、确定节点间断有限元方法，根据步骤s1得到的目标连续变化介质模型参数确定步骤s2得到的波动方程的空间离散格式；

7、s4、确定时间离散方法和时间步长计算方法；

8、s5、将步骤s3得到的波动方程的空间离散格式代入步骤s4得到的时间离散方法中，并根据步骤s4得到的时间步长进行迭代，计算测试时长内每个时间步长下的解，构成最终的波场正演结果。

9、优选地，步骤s1中，网格信息包含单元节点坐标和单元位置关系。

10、优选地，步骤s2中，波动方程为一阶形式的声波方程和弹性波方程，具体如下：

11、

12、其中，u为解向量，g为震源项，矩阵a1和a2为系数矩阵，t为时间，x，y为二维空间坐标。

13、优选地，步骤s3具体为：

14、s301、将计算域划分为互不重叠的三角形单元后，利用插值多项式对单元内局部解进行近似；

15、s302、将测试函数与波动方程相乘并将结果在单元内积分，获得待求解的方程系统，通过分部积分法获得方程系统的弱形式；

16、s303、定义参考单元dr，对参数矩阵进行拉格朗日多项式近似，并进行坐标变换，得到参考坐标系下方程系统的弱形式，将边界积分代入方程系统的弱形式得到最终的离散格式。

17、更优选地，步骤s301中，每个单元内的局部解u(x,t)为：

18、

19、其中，li(x)是基于单元dk内部插值锚点的拉格朗日插值多项式，np为插值点个数。

20、更优选地，步骤s302中，方程系统的弱形式如下：

21、

22、其中，f为边界处的流通量，u为解向量，g为震源项，矩阵a1和a2为系数矩阵，t为时间，x，y为二维空间坐标，li为测试函数，dk为第k个三角形单元所包含的积分区域，为三角形单元的边界。

23、更优选地，步骤s303中，离散格式为：

24、

25、其中，uj为单元内锚点处的解向量，jk为单元dk的坐标变换雅可比行列式，li为测试函数，lj为单元内局部解的插值多项式，aξ,m为，ζ＝(ξ,η)为参考坐标系，lm为系数矩阵的插值多项式，gj为震源向量值，为单元边界处的雅可比行列式，te为坐标变换矩阵，与为系数矩阵，为单元边界锚点处的解向量，dr,e为参考单元的第e条边界边界。

26、优选地，步骤s4中，时间离散方法采用四阶runge-kutta方法，方程y'＝f(t,y),y(t0)＝y0有如下形式：

27、

28、其中，yn为tn时刻下方程的解，k1，k2，k3，k4为tn与tn+1间各个节点处函数f(t,y)的值，δt为时间步长。

29、优选地，步骤s5中，数值误差e具体为：

30、

31、其中，n为计算结果总时间步数，sndg为本发明方法计算结果，ssem为被广泛认可的sem方法计算结果。

32、第二方面，本发明实施例提供了一种基于节点间断有限元方法的波场正演模拟系统，包括：

33、参数模块，对目标连续变化介质模型进行非结构网格划分，导出网格信息并关联相关物性参数，得到目标连续变化介质模型参数；

34、方程模块，确定波动方程；

35、离散模块，确定节点间断有限元方法，根据参数模块得到的目标连续变化介质模型参数确定方程模块得到的波动方程的空间离散格式；

36、计算模块，确定时间离散方法和时间步长计算方法；

37、模拟模块，将离散模块得到的波动方程的空间离散格式代入计算模块得到的时间离散方法中，根据得到的数值误差分析方法精确性。

38、与现有技术相比，本发明至少具有以下有益效果：

39、一种基于节点间断有限元方法的波场正演模拟方法，分为如下几个步骤：建立模型与网格划分；确定波动方程；确定空间离散方法；确定时间离散方法；数值求解以及方法验证。方法步骤清晰明确，各步骤间关系紧密，逻辑清晰，便于实现。

40、进一步的，步骤s2中，网格信息包含单元节点坐标，其目的在于可以通过单元节点坐标计算得到单元内锚点的具体坐标，根据锚点坐标和连续变化介质模型可以确定锚点处的物性参数，为后续系数矩阵的多项式近似提供依据。网格信息包含单元位置关系，其目的在于可以通过单元间的位置关系在后续系统实现中恢复网格。

41、进一步的，步骤s301利用插值多项式对互不重叠的三角形单元内的局部解进行近似，其目的在于得出局部解的离散格式，用于步骤s302中方程弱形式的推导。步骤s303中对参数矩阵进行多项式近似，其目的在于，将参数矩阵以离散形式表示，使得矩阵能够提出积分号外，进一步的将锚点处的物性参数加入计算，使方法能够处理连续变化介质。步骤s303中定义参考单元及参考坐标系，并通过坐标变换将方程系统由非参考三角单元转换至参考单元中统一计算，便于计算机并行计算，提高计算效率。

42、进一步的，步骤s4中采用四阶runge-kutta方法作为时间离散方法，其优势在于，该方法间接采用了taylor展开方法，使用数个节点处的函数值的线性组合替代了函数的导数，避免了求导计算的同时保证了方法的高精度，形式简单便于实现。

43、可以理解的是，上述第二方面的有益效果可以参见上述第一方面中的相关描述，在此不再赘述。

44、综上所述，本发明模拟方法计算效率高，所得的正演结果精度高。

45、下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。