显格式有限差分具有算法实现简单、内存需求低和计算量小等性能,因而在求解地震波场数值模拟领域中应用十分广泛。然而,其最大时间步长严格受制于Courant-Friedrichs-Lewy (CFL)稳定性条件限制,这要求我们在小尺度结构和高速模型中必须采用更小的网格才能获得稳定的模拟结果,因此其计算量势必会大幅增加。中科院地质与地球物理所地球与行星物理院重点实验室的博士后高英杰与张金海研究员、合作导师姚振兴研究员,采用特征值扰动法,有效破除了显格式有限差分的CFL稳定性条件,使得我们可以采用比以往大得多的时间步长进行稳定的数值模拟。
当时间步长小于CFL约束的最大时间步长时,迭代矩阵的特征值的模始终位于单位圆上或者单位圆内,如图1(a)和(b)所示;当时间步长超过CFL约束的最大时间步长时,迭代矩阵特征值的模会有一部分超出单位圆,从而导致迭代失稳,如图1(c)所示。特征值扰动法可以将不稳定特征值的模归一化到单位圆上,如图1(d)所示,使其仍然满足特征值不大于1的基本条件,从而可以利用任意时间步长实现稳定的数值求解。然而,采用大步长模拟将会面临严重的时间数值频散,需要采用逆时间离散变换方法才能获得精确的模拟结果,如图2(b)所示。最大时间步长可以是稳定性条件所允许步长的4~5倍,并且达到了地震波高截频的奈奎斯特采样上限,从而为深部探测和长时程探测提供了精度不受影响但计算效率更高的新方法,该技术对于全球尺度的地震波模拟也具有重要意义。
图1 特征值扰动示意图。采用不同时间步长所获得的特征值:(a) Δt = 0.5 ms;(b) Δt = 1.5 ms;(c) Δt = 2 ms;(d) Δt = 2 ms。(d)为(c)经过特征值扰动后的结果。本实验中,CFL所允许的最大时间步长为1.53 ms
图2 利用不同时间步长获得的单道波形记录。依次测试了时间步长Δt = 2, 3, 4, 5, 6, 7 ms。 (a)是采用过特征值扰动算法之后得到的模拟波形;(b)是对(a)采用了逆时离散时间变换后的结果。本实验中,CFL所允许的最大时间步长为1.53 ms,震源子波函数对应的奈奎斯特采样上限为6.7 ms
研究成果以Letter的形式发表于Geophysics。(Gao Y, Zhang J, Yao Z. Removing the stability limit of the explicit finite-difference scheme with eigenvalue perturbation[J]. Geophysics, 2018, 83(6): A93–A98.)(原文链接)