地震学成像方法是地球内部结构探测的重要手段。地震波场模拟是逆时偏移、全波形反演等地震波成像方法的基础，高精度、高效率的地震波场模拟方法直接影响成像的质量和效率。近年来，谱元法因其高精度和高效率的特点在地震学中得到了广泛的应用。谱元法通过在每个单元上采用高阶多项式展开，也被称为高阶的有限元法或区域分解的谱方法，它兼顾了有限元法的灵活性和谱方法的指数收敛性。由于重合的插值节点与积分节点可以获得对角的质量矩阵，避免了求解大型线性代数方程组，因而具有较高的计算效率。然而，对角的质量矩阵依赖Gauss-Lobatto（GL）积分。 

　　目前，谱元法多应用于四边形、六面体单元。由于这些单元自身的正交性质，高维单元的GL积分通过一维张量积实现。尽管四边形网格谱元法具有较高的计算精度，然而它在刻画几何模型方面远没有三角形单元灵活。但是三角形单元为非正交单元，不能采用一维GL张量积对波动方程进行离散。传统的三角网格谱元法采用广义Newton-Cotes积分对波动方程进行离散，该积分仅具有N阶积分精度，远不能近似2N阶的质量矩阵和2N-2阶的刚度矩阵，使得传统的三角网格谱元法不具有指数收敛性。 

　　为此，中国科学院地质与地球物理研究所刘有山博士后和徐涛研究员等人发展直接应用于三角网格中的GL积分方法。发展一种基于p范数的目标函数，获取了三角形中的高阶（7-9）阶GL积分公式，该求积公式适用于任意的二阶偏微分方程，我们将其初步应用于三角网格谱元法的求解中。数值算例表明：在计算精度方面，基于优化求积点的三角形网格谱元法将传统的三角网格谱元法的计算精度提高近一个数量级，使之能与高精度的四边形网格谱元法相当；在计算效率方面，基于优化求积点的三角形网格谱元法远高于传统的三角网格谱元法，略低于四边形网格谱元法。 

图1 t=0.75 s时刻的垂向位移波场快照
(a) 10阶传统的三角网格谱元法(TSEM)；(b) 14阶传统的三角网格谱元法；
(b) (c) 7阶四边形网格谱元法(QSEM)；(d) 7阶优化点积的三角网格谱元法(OTSEM) 　　  

　　上述成果近期发表在国际期刊JCP上（Liu et al., 2017. Higher-order triangular spectral element method with optimized cubature points for seismic wavefield modeling, Journal of Computational Physics, 336: 458-480）。该研究得到科技部国家重点研发计划资助（2016YFC0600101）。