本发明是一种应用于复杂地表、复杂构造中地震波波场数值模拟技术，具体涉及一种有限元类数值求解算法。
背景技术：
目前，在我国南方和西部地区，油气地震勘探的重点正转向丘陵、山前构造带等复杂地区。这些地区地表条件异常复杂，地形起伏剧烈，高差变化非常大，岩性速度变化大导致近地表结构不均匀性严重，同时地下构造复杂，如褶皱强烈、断层发育、构造陡峻、地层变化大等。它们导致了这些地区地震资料信噪比低及静校正难等各方面的问题。从根本上解决这些勘探问题，需要我们对起伏地表条件下地震波的传播规律和波场特征进行深入的认识，而有限元法是进行复杂地表和复杂构造地区地震波数值模拟最有效的技术手段。
间断galerkin有限元法(dgfem)是一种高阶有限元法。基于数值流理论的dg-fem本质上是有限元法和有限体法的结合，在单元内部使用有限元法处理，单元边界上采用了有限体法中数值流的处理思想，因此它继承了有限元法的高阶和有限体法的局部特性等优点，同时也克服了这两个方法的一些缺陷，能够实现高精度、低频散、有效的数值模拟。它可以使用非结构网格单元(包括三角形或四面体网格)，能够根据介质的分布特征设计出最优的网格，因此，间断galerkin有限元法特别适应起伏地表及复杂构造条件下地震波传播数值模拟。
dg-fem由于本身高阶收敛等特性，当采用一般的显式时间格式时，对时间步长的要求比较严格，例如在龙哥库塔时间格式中，随着空间阶数p的增加，步长以1/(2p+1)的比例降低。另外，当使用有限元网格对地下介质模型进行精细剖分时，不可避免的出现“小网格高速体”单元，根据cfl稳定性条件，对步长的限制更加突出。当模型中存在局部较密网格时，虽然绝大部分区域网格单元较大，当采用全局时间步长时，为了满足局部密网格严格的时间步长，全局时间步长会很小，这样会严重制约整个模拟计算的效率。
克服显式时间格式步长限制的通常一个方法是采用隐式时间格式，但一般隐式时间格式需要对大规模稀疏矩阵求逆，这并不是一件容易的事。对于实际中使用dg-fem模拟弹性波传播，即使采用稀疏存储，对离散后的大型稀疏矩阵的存储也是难以承受的，三维情况下更加突出。另一个可取的方法是采用局部时间步长，但其完成比较复杂。
相比高阶的间断galerkin有限元法，低阶方法允许采用更大的时间步长。间断galerkin有限元法拥有良好的局部特性，它的求解依赖于逐个单元，它允许每个单元选择独立的空间阶数进行计算，因此可以考虑对拥有更严格时间步长限制的单元采用低阶模拟，而在其它单元采用高阶模拟保证模拟精度，即所谓的阶数自适应或p自适应。
技术实现要素：
本发明针对使用间断galerkin有限元法(lf-dg)进行地震波传播数值模拟时，模型中局部小网格单元或畸形单元会导致较小全局时间步长的问题，通过利用间断galerkin有限元法的局部特性对p自适应算法进行实现，允许各个单元设置独立的空间阶数，进而获得优化的局部时间步长，提高模拟计算效率。
本发明的目的及解决其技术问题是采用以下技术方案来实现的。
所述方法包括以下步骤：
(1)定义网格剖分密度，按网格密度对计算区域进行非结构网格剖分；
(2)确定每个网格单元上的多项式展开阶数p；
(3)对相邻的两个网格单元解耦，选择合适的数值流通量格式作为相邻单元间波场交换方式；
(4)对解耦方程数值化，获得每个网格单元上的空间离散方程；
(5)求相邻两个单元公共边界上波场函数之间的内积，获得边值映射矩阵；
(6)选择合适的时间积分格式；
(7)利用边值映射矩阵对不同阶数的相邻单元波场进行映射，计算空间离散部分，施加震源条件，更新波场分量，获得模拟结果。
其中步骤(2)是通过以下步骤获得：
(1)输入模拟时间步长dt和最大模拟阶数pmax；
(2)确定每阶算法的稳定性条件数cfl(p)；
(3)计算每个单元的时间步长标准；
(4)计算pmax下，每个单元k的稳定性条件；
(5)对dtk≥dt所有单元采用pmax阶；
(6)依次重复(4)(5)的步骤，依次筛选出阶数pmax-1、pmax-2…下的计算单元；
(7)当所有单元没有获得确定的计算阶数时，提示减小模拟时间步长大小dt。
步骤(5)是通过以下步骤得出：
(1)将相邻两单元格的函数映射到单元边界上，设它们在边界上分别由两组不同阶数的正交模态空间
(2)推导不同阶数展开的函数
(3)用函数在单元上的模态空间坐标及边界积分算子表示函数在单元边界上的内积
所述时间积分格式为龙哥库塔时间格式、蛙跳时间格式。
所述数值流为中心数值流、局部lax-friedrich。
步骤(2)中imn表示维度为(pm+1)*(pn+1)的单位矩阵。
其中
函数在边界上的物理空间坐标可以由函数在单元上的模态空间坐标表示为
其中边界内积算子
当所有单元没有获得确定的计算阶数时，提示减小模拟时间步长大小dt。
本发明有以下优点和有益效果：本发明适用复杂条件(特别是剧烈地形起伏等)中地震波传播的高精度高效数值模拟，通过利用间断galerkin有限元法的局部特性实现了p自适应算法，允许各个单元设置独立的空间阶数，进而获得优化的局部时间步长，提高模拟计算效率。该发明理论新颖，技术流程实用性强，适用于业界广泛运用。
附图说明
在下文中将基于实施例并参考附图来对本发明进行更详细的描述。其中：
图1是均匀模型示意图及其网格分布示意图
图2(a)是空间阶数分布图
图2(b)是采用时间步长为dt＝0.20ms时，p自适应模拟结果效果图
图2(c)是采用时间步长为dt＝0.20ms时，单阶模拟结果差异效果图
图3(a)是空间阶数分布图
图3(b)是采用时间步长为dt＝0.25ms时，p自适应模拟结果效果图
图3(c)是采用时间步长为dt＝0.25ms时，单阶模拟结果差异效果图
图4(a)是空间阶数分布图
图4(b)是采用时间步长为dt＝0.40ms时，p自适应模拟结果效果图
图4(c)是采用时间步长为dt＝0.40ms时，单阶模拟结果差异效果图
图5是l2误差及计算效率对照表
在附图中，相同的部件使用相同的附图标记。附图并未按照实际的比例。
具体实施方式
下面将结合附图对本发明作进一步说明。
本发明本发明针对使用间断galerkin有限元法进行地震波传播数值模拟时，模型中局部小网格单元或畸形单元会导致较小全局时间步长的问题，通过重新考查两个不同阶数的波场函数在同一积分区域上(单元边界)的内积，充分应用配点微分法求内积的方式，巧妙推导了不同阶数展开波场之间的映射矩阵，实现了p自适应模拟，并提供了为各个单元设置独立的空间多项式展开阶数的方法，进而获得优化的局部时间步长，提高模拟计算效率。该发明的具体步骤包括：
根据模型速度分布以及地震波有限元法的空间采样需求定义网格密度函数
输入模拟计算时间步长及采用的最大空间阶数pmax，根据模拟的稳定性准则设置每个单元上最合适的空间多项式展开阶数。确定网格单元的空间多项式展开阶数的方法步骤见附图1。
应用间断galerkin有限元法的基本原理推导弹性波一阶速度-应力方程对应的积分方程(有限元方程)，实现相邻网格单元的解耦，其中选择合适的数值流作为相邻单元间的波场交换方式，如中心数值流、局部lax-friedrich等。
采用配点微分法求内积的方式将解耦积分方程数值化，获得每个单元上的空间离散方程。p自适应算法允许各个单元采用独立的空间多项式展开阶数。只有当求波场函数在单元边界上内积项时，相邻单元采用的空间阶数才会影响本单元波场的求取，而其他体积分项的求取只与本单元的空间阶数选取有关，因此重点考察相邻两个单元公共边界上波场函数的内积，表示为
其中，
推导不同阶数展开的函数
其中，imn表示维度为(pm+1)*(pn+1)的单位矩阵。由边界上模态空间和物理空间之间的映射关系，有
另外从整个单元上看，由函数在边界上的物理空间坐标可以由函数在单元上的模态空间坐标表示为
综合上面的两组表达式有
最终，函数在单元边界上的内积可以由函数在单元上的模态空间坐标及边界积分算子表示为(这里直接采用<2-40>式中的节点形式表示)
其中，边界内积算子
选择合适的时间积分格式，如龙哥库塔时间格式、蛙跳时间格式等。
利用映射矩阵tmn对不同阶数的相邻单元波场进行映射，计算空间离散部分，施加震源条件，更新波场分量，获得最终模拟结果。
实施例：
附图2为均匀模型示意图及网格分布式情况；
模型采用p3阶算法模拟。为了满足全局网格对稳定性的要求，时间步长dt≤0.11ms，采用0.1ms的模拟时间步长。当采用dt＝0.2ms时间步长模拟时，从附图3中可以看到为了满足稳定性要求，局部网格加密处有少部分单元需要采用p2阶单元计算，其他区域仍采用p3阶单元模拟，附图3(b)为p自适应模拟结果，附图3(c)为p自适应模拟结果与单一p3阶模拟结果的差异。比较数值量级，模拟结果的峰值达到4*104，差异的峰值为1000，即差异约为2-5％。
附图4中当采用dt＝0.25ms时，局部网格加密处更多单元需要采用p2阶模拟。与单一p3阶模拟结果差异约为3-7.5％。
附图5中当采用dt＝0.4ms时，局部网格加密处大部分单元需要采用p1-p2阶模拟，模型其他区域也有部分单元需要采用p2阶模拟。与单一p3阶模拟结果差异约为7.5％。
从图5中，可以看到p自适应算法能够在满足精度要求的同时，通过增加模拟时间步长，明显提高计算效率。
虽然已经参考优选实施例对本发明进行了描述，但在不脱离本发明的范围的情况下，可以对其进行各种改进并且可以用等效物替换其中的部件。尤其是，只要不存在结构冲突，各个实施例中所提到的各项技术特征均可以任意方式组合起来。本发明并不局限于文中公开的特定实施例，而是包括落入权利要求的范围内的所有技术方案。