地震波是探测地球内部结构的重要手段之一，基于走时数据的全球层析成像的结果已经获得了全球尺度的认识，其中最为显著的地震波速度结构是在非洲和太平洋板块下方的两个大尺度低速区。随着高性能计算机集群和并行算法的发展，复杂介质地震波谱元法数值模拟快速发展。新一代基于地震波波形的较高分辨率的全球波形层析成像得到了更加清晰、趋于统一的大尺度地幔结构（几百公里到千公里）。虽然更高分辨率（数公里到百公里）的小尺度结构体的信息隐藏于全球波形数据中，但受制于全球波形正反演计算能力，高频信息不得不被舍弃，从而未能由全球层析成像反演出来。从这点上说，计算方法落后于成像理论的发展。 地球深部关键区域波形层析成像（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）。
图3显示了如果目标模型是光滑且各向同性的，在不同的初始模型下，反演所得的模型均收敛于光滑的目标模型。图4显示如果目标模型是多尺度且各向同性的，基于有限带宽的波形，各向同性反演（图4d）所得模型无法收敛于真实的光滑目标模型（图4b）。而各向异性反演所得的速度和各向异性模型（图4e和图4f）也同样无法收敛于光滑目标模型（图4b和图4c）。但图5显示对于震源和台站均在反演区域外的路径而言（SA），在图4光滑目标模型bc和各向异性反演所得模型ef中计算所得波形完全一致；而对于震源在外侧，而台站在反演区域里的路径而言（SB），图4b、图4c和图4e、图4f两模型中计算所得波形却不一致，说明反演所得模型整体发生了坐标变换(particle relabeling transformations，David Al-Attar et al., 2016)。
图3 光滑各向同性声波介质目标模型下，不同初始模型所得的全波形反演结果。（a）光滑的目标模型以及震源（红色五角星）和台站（黑色三角）分布。（b）目标模型的一种特殊光滑（homogenization，Capdeville et al., 2018）。（c、f）两种不同的初始速度模型。（d、g）两种不同初始模型下全波形反演所得速度模型。（e、h）两种全波形反演模型（d和g）对应的光滑模型，和光滑的目标模型（b）一致
图4 多尺度目标声波介质模型下，全波形反演的结果。(a)多尺度目标模型。(b)目标模型的一种光滑模型。(c)多尺度目标模型光滑以后所得到的等效各项异性。(d)各项同性全波形反演所得速度模型。(e)各向异性全波形反演所得速度模型。(f)各项异性全波形反演所得各向异性
图5 图4中SA和SB两条路径下的波形对比图。红线是反演所得模型和初始模型下数值模拟波形之差，蓝线则是是目标模型和初始模型下数值模拟波形之差。（SA）震源和台站均在反演区域外；（SB）震源在外侧，而台站在反演区域里
上述地震波的混合数值模拟，超高阶谱元法的正演计算，以及非均一化全波形反演的非唯一性的研究，为下一代子区域全波形层析成像奠定了重要基础，可用于探索位于地球深部任意位置的多尺度地震波速度结构。
研究成果分别发表于国际学术期刊Journal of Geophysical Research: Solid Earth, Geophysics, Geophysical Journal International。该研究得到国家自然科学基金（41625016，41888101，4200404）项目，国家留学基金委（201804910289）项目和法国国家研究署（ANR-16-CE31-0022-01）的共同资助。
1. 吕超, 赵亮*, Capdeville Yann. Novel hybrid numerical simulation of the wave equation by combining physical and numerical representation theorems and a review of hybrid methodologies [J]. Journal of Geophysical Research: Solid Earth, 2022, 127: e2021JB022368. DOI: 10.1029/2021JB022368.
2. 吕超, Yann Capdeville, 吕刚, 赵亮*. Removing the Courant-Friedrichs-Lewy stability criterion of the explicit time-domain very high degree spectral-element method with eigenvalue perturbation [J]. Geophysics, 2021, 86: T411-T419. DOI: 10.1190/geo2020-0623.1.
3. 吕超, Yann Capdeville*, David Al-Attar, 赵亮. Intrinsic non-uniqueness of the acoustic full waveform inverse problem [J]. Geophysical Journal International, 2021, 226(2): 795–802. DOI: 10.1093/gji/ggab134.
