黏弹性波动方程求解的四阶时间离散立体建模方法

《Petroleum Science》：Stereo-modeling methods with fourth-order time discretizations for solving the viscoelastic wave equation

　　

地下介质普遍表现出黏弹性特征。基于黏弹性波方程开发正演模拟算法对研究波传播特性至关重要。正演模拟结果的精度取决于时间与空间离散化。常规有限差分(CFD)方法在时间上通常达到二阶精度，但在复杂黏弹性波模拟中难以兼顾精度与计算效率，其中S波更易受数值频散影响。因此，实际应用中亟需稳定、高阶时间离散方法结合高精度空间离散。本文针对达朗贝尔(D’Alembert)波方程，在不同理论框架下发展了两种立体建模(STEM)方法，实现时间与平面均四阶精度。具体而言，空间离散采用基于方向导数的四阶STEM算子。研究人员进而结合STEM算子发展两种四阶时间精度方法：牛顿体系下的近解析中心差分法(NACD)与伯克霍夫(Birkhoffian)体系下的辛立体建模法(SSM)。通过数值频散、效率、误差、能量稳定性及正演性能等一系列实验分析NACD与SSM的数值精度与长时间传播稳定性。精度方面，所提四阶时间方法优于传统二阶方法；即便与四阶时间的Lax-Wendroff修正法(LWC)相比，NACD与SSM更贴合解析解，更准确捕捉黏弹性介质衰减，以更高效率更好抑制数值频散。NACD与SSM精度相当，但NACD计算效率更高。长时间稳定性方面，相比同算子的二阶NADM，NACD与SSM长期误差积累更小；稳定性测试表明SSM因辛时间方法更好保持能量稳定。故NACD更适常规地震勘探，SSM更适长时间模拟如大规模模拟与多次波成像。

研究背景方面，地下储层普遍具黏弹性与衰减特征，基于全波方程的正演是合成地震记录与研究波传播的关键步骤。常规有限差分(CFD)方法时间上多为二阶精度，粗网格或大速度反差下空间高阶算子会带来数值频散、边界处理困难与内存膨胀；若加密网格则时间步需缩小，计算成本指数上升。时间离散研究相对滞后，典型二阶时间离散在常用步长下频散误差显著。已有四阶时间Lax-Wendroff修正法(LWC)通过方程将时间导数转为空间导数以提升精度，而立体建模(STEM)算子借位移与梯度联合截断泰勒展开抑制频散且紧致。牛顿体系下时间推进不考虑能量守恒，长时间模拟易有人工耗散；辛方法在哈密顿(Hamiltonian)体系可保结构，但黏弹性方程须纳入更广义的伯克霍夫(Birkhoffian)体系。针对上述矛盾，研究人员在《Petroleum Science》发表本文，目标是在达朗贝尔型黏阻尼波方程上，耦合四阶STEM空间算子与四阶时间离散，分别给出牛顿体系NACD与伯克霍夫体系SSM，并系统比较精度与长时间稳定性，为复杂介质正演与成像提供基础。

主要关键技术方法上，空间离散采用基于方向导数的四阶STEM算子，同时利用位移u₁、u₃及其梯度构造高阶偏导逼近，紧致三点支撑抑制各向异性频散。时间离散分两条线：其一在牛顿体系对u₁分量作时间导数递推，保留至四阶时间项并结合中心差分得NACD；其二将黏弹性方程改写为广义Birkhoffian形式，采用Qin和Wang推导的四阶辛分割龙格库塔(PRK)系数构造含四中间阶段的SSM，并证其辛性。稳定性分析借谐波代入求增长矩阵特征值，用遍历法估Courant数上限。精度与频散评测以解析解（Aki等）为基准，定义频散比R=v/v_num及相对误差E_r，并在两层强界面模型与均匀长时模型上做快照、波形与相对能量H(u)演化对照。

研究结果部分，各小节结论如下。稳定性准则小节：经谐波展开得到NACD、SSM与LWC增长矩阵，数值遍历得稳定条件η≤0.55(NACD)、η≤0.42(SSM)、η≤0.50(LWC)；SSM因单步内四阶段重复作用空间算子，对高波数更敏感故稳定性更严，NACD形式紧凑限制较松。解析解对比小节：12×12 km²粗网格(h=0.04 km, Δt=4 ms)下，CFD与LWC快照可见明显频散，NACD与SSM的u₁、u₃分量波形与解析解吻合最好，验证四阶时空精度。数值频散分析小节：取η=0.2、S_p扫描，CFD最大频散误差27.0%，LWC 24.3%，NACD约10.8%，SSM约10.6%；S_p=0.4极坐标下NACD最大误差0.0506、SSM 0.0416，SSM各向异性稍优，二者均压制S波频散。长时间传播稳定性小节：201×201网格模拟60 s，二阶NADM误差持续累积发散，NACD与SSM稳定，末态SSM误差比NACD低约10%；1D相对能量H(u)演化中SSM仅微幅波动保能量，NACD因泰勒展开时间项出现明显能量起伏，证SSM长时间优。

正演建模小节：10×10 km²两层强界面模型(h=0.025 km)中，LWC粗网格低频层与界面处频散严重，NACD与SSM几乎无频散且清晰分辨直P/S、反射PP/PS/SS/SP、透射Pp/Ps/Ss/Sp等转换波；将LWC网格加密至h=0.010 km可达相近画质但CPU 1364.91 s，NACD仅135.31 s、SSM 630.95 s，存储分别为LWC的48.14%与64.19%；单点FLOPs NACD 462、SSM 892高于LWC 290，但因总网格数骤减，总FLOPs NACD 0.6×10¹¹、SSM 1.1×10¹¹远小于LWC 5.8×10¹¹。黏弹与纯弹性对照显示NACD/SSM正确再现振幅衰减与频谱衰减。

讨论与结论翻译：研究人员指出，NACD与SSM在四阶STEM空间算子上均达四阶时空精度，比二阶CFD与四阶LWC在粗网格更准、更高效、频散更低；NACD计算量最小适常规勘探，SSM因Birkhoffian辛结构长时间能量稳定适大尺度与多次波成像；D’Alembert方程非严格常Q模型，极低速低Q宽频仍需改进。该研究意义在于为黏弹性介质正演提供两套可选高阶方案，并示范牛顿与Birkhoffian双体系耦合STEM算子的建模范式，支撑后续逆时偏移与黏弹性全波形反演。