地震波是探测地球内部结构的重要手段之一，基于走时数据的全球层析成像的结果已经获得了全球尺度的认识，其中最为显著的地震波速度结构是在非洲和太平洋板块下方的两个大尺度低速区。随着高性能计算机集群和并行算法的发展，复杂介质地震波谱元法数值模拟快速发展。新一代基于地震波波形的较高分辨率的全球波形层析成像得到了更加清晰、趋于统一的大尺度地幔结构（几百公里到千公里）。虽然更高分辨率（数公里到百公里）的小尺度结构体的信息隐藏于全球波形数据中，但受制于全球波形正反演计算能力，高频信息不得不被舍弃，从而未能由全球层析成像反演出来。从这点上说，计算方法落后于成像理论的发展。
地球深部关键区域波形层析成像（Box tomography），结合了高效的地震波混合数值模拟方法和全波形反演理论（Full waveform inversion）的优点,成像所需的大量正反演计算只在子区域内部进行，可以用较少的计算资源对地球深部关键区域进行高频全波形反演成像，因此Box tomography被认为是应用全波形反演方法寻找小尺度结构体的下一代技术。而作为其正演基础的“地震波混合数值模拟”的计算效率则是关键所在。
针对地震波混合数值模拟方法，中国科学院地质与地球物理研究所岩石圈演化国家重点实验室吕超博士后与合作导师赵亮研究员，法国科学院Yann Capdeville教授展开合作，对现有的5种地震波混合数值模拟方法进行理论和数值分析对比评述，通过将地震波经典的物理形式“表示定理”（Representation theory）和其等效的数值表达形式相结合，提出了一种全新、灵活和高效的混合数值模拟技术。
该技术主要分为两步：第一,在全球地震波正演计算过程中只需要存储下子区域边界处（图1a中的绿线）的位移、应变和加速度三个物理量；第二,将存储的三个物理量转换成等效体力加载到子区域的同样的边界处进行混合地震数值模拟。
该新方法第一步中的物理量的计算存储，第二步中等效体力的加载均只与子区域的边界有关，混合计算所需的存储量有降维优势，且所得混合波形精度非常精确，可通过调整空间单元数量及单元内部的阶数来控制（见图1b中的波形对比）（Lyu et al., JGR, 2022）。
图1 混合数值模拟波场快照及波形图。图a显示了在混合数值模拟计算中，将获取的物理量加载进子区域中的过程。因子区域内部存在小尺度异常体（灰色的圆盘），只有残差波场传出到子区域外部。黑色五角星代表远端的震源，黑色倒三角代表子区域内部的台站。图b为波形对比，黑线代表小尺度异常体模型下，全局正演计算所得的波形；红虚线代表包含小尺度异常体子区域模型下，混合计算所得波形；绿线所示10000倍放大了黑线和红虚线波形差；蓝虚线代表没有小尺度异常体的全局正演计算所得波形
就正演计算而言，常用的谱元法数值模拟所采用的网格单元单个方向上往往只采用5个Gauss-Lobatto-Legendre（GLL）格点，主要有两个原因：1）所有物性间断面必须采用显示的网格划分以确保数值模拟精度；2）超高阶谱元法的最小网格间距太小使得可用时间步长太小，无法进行长时程数值模拟。为尝试解决此问题，研究团队结合了正向时间离散变换(FTDT)、本征值扰动法(Perturbation)和反向时间离散变换(ITDT)，首次突破了显式时间域超高阶谱元法的稳定性条件，图2显示混合前后波场快照对比图。本征值扰动法通过剔除不稳定特征值突破了CFL约束，使得可用的时间步长接近奈奎斯特采样极限，极大地减少了正演计算所需时间迭代总次数。虽然大时间步长模拟将会产生严重的时间数值频散，正反时间离散变换的结合使用则确保了数值模拟的精度（Lyu et al., Geophysics, 2021）。
图2 波场快照对比图。(a)超高阶谱元法数值模拟在4个不同时间步长下的波动方程数值模拟。(b)结合了正向时间离散变换(FTDT)、本征值扰动法(Perturbation)和反向时间离散变换以后的超高阶谱元法数值模拟在4个不同时间步长下的波动方程数值模拟
针对多尺度全波形反演的非唯一性问题而言，基于有限频带波形数据，多尺度弹性波介质（固态地球）的反问题被证明具有唯一性，即全波形反演能收敛于目标模型（Capdeville et al., 2018）。为研究多尺度声波介质的全波形反演（海洋及液态外核）的非唯一性，研究团队首次在数值上证明了光滑的各向同性目标模型的反问题具有唯一性（图3），而多尺度声波介质的反问题却具有内在非唯一性（图4和图5），且非唯一性可能源于一种坐标变换—质点重标记变换(particle relabeling transformations, David Al-Attar et al., 2016) （Lyu et al., GJI, 2021）。
