SeisFlows--Flexible waveform inversion software

SeiFlows————灵活的波形反演软件

摘要

SeisFlows 是一个开源 Python 包，为石油和天然气勘探、地震层析成像、医学成像和其他领域的研究提供可定制的波形反演工作流程和框架。通过继承默认的反演或偏移类，新方法可以在 SeisFlows 中快速原型化，并且可以在应用到更昂贵的 3D 问题之前在 2D 示例上测试代码。必须使用外部软件包（例如 SPECFEM3D）执行波形模拟。与外部解算器接口的能力提供了灵活性，并且选择 SPECFEM3D 作为默认选项提供可选的 GPU 加速和其他有用的功能。通过支持高性能计算 (HPC) 系统的大规模并行求解器和接口，可以执行具有数千条地震轨迹和数十亿模型参数的反演。到目前为止，SeisFlows 已在由国防部、雪佛龙公司、道达尔公司、普林斯顿大学和阿拉斯加大学费尔班克斯分校管理的集群上运行。

[TOC]

1. 介绍

波形反演是一种强大且计算量大的方法，用于估计物体的材料属性和成像其内部结构（Virieux 和 Optero，2009；Burstedde 和 Ghattas，2009）。 无论是应用于人体、油气藏还是在大陆克拉通，需要仔细集成软件组件，包括波动方程求解器、非线性优化库和信号处理程序。

许多波形反演包是石油和天然气公司开发的专有代码。 通常，独立研究人员无法使用此类软件。

在工业之外，已经开发了许多开源波形反演包。 MADAGASCAR 是一个广泛使用的勘探地球物理软件包，提供了多种功能，包括一些波形反演功能（Fomel 等人，2012 年）。 SEISCOPE 是另一种广泛使用的工具，它提供 FORTRAN 90 非c线性优化例程，而不是完整的反演工作流程（M_etivier 和 Brossier，2016 年）。 最近的三个软件包 PySIT、LASIF 和 SimPEG 使用 Python 进行软件集成和数据处理任务。 PySIT 以工具箱格式提供强大的波形反演功能（Hewett 和 Demanet，2013 年）。 LASIF 为组织和可视化地震断层扫描结果提供了一个有用的系统（Krischer 等人，2015a）。 SimPEG 为地球物理反演提供了广泛的框架，并附带基于有限体积法的模拟工具（Cockett 等人，2015 年）。

SeisFlows 与以前的软件不同，它提供了一个强调灵活性和 HPC 可移植性的高级反演工作流程。 为了在广泛的社区中使用并在不断发展的研究领域中保持相关性，软件必须在硬件和软件环境中具有灵活性和可移植性。 对于方法论研究，希望能够在工作站或小型集群上进行原型设计。 对于大型 3D 应用程序，在大规模并行环境中运行是必不可少的。 由于集群之间缺乏标准化，可移植性在大规模应用中尤其困难

为了帮助实现这些设计目标，SeisFlows 将反演问题抽象为六个部分：(1) 求解器、(2) 预处理、(3) 后处理、(4) 非线性优化、(5) 系统和 (6) 工作流。 SeisFlows 底层的源代码基于这些类别以模块化方式构建，并且为每个用户提供了多种选择。 例如，如果地震断层扫描项目中的研究区域扩大，用户可以将笛卡尔求解器换成球形全地球求解器。 如果 PBS 集群下线，而 SLURM 集群上线替换它，用户可以将 PBS 系统接口换成 SLURM 系统接口。 如果主包中缺少所需的功能，用户可以贡献自己的类或重载默认类。

名称	内容
工作流	反演 成像
系统	串行 多线程 批处理 分布式集群管理系统 大型计算节点集群的高度可伸缩和容错的集群管理器和作业调度系统
求解器	SPECFEM2D SPECFEM3D SPECFEM3D_GLOBE
非线性优化	梯度下降 非线性共轭梯度法 拟牛顿算法
预处理	波前 包络 瞬时相
后处理	通过平滑正则化 Tikhonov 总变换

图 1. SeisFlows 包的结构。 为组成该包的六个部分中的每一个列出了一系列主要选项或设置。

本文分为四个部分。 首先，我们详细阐述了 SeisFlows 背后的动机。 其次，简要回顾了波形反演的理论。 第三，我们描述了六个组成部分中的每一个，以及它们如何协同工作来管理反演的复杂性。 最后，我们通过一个勘探地球物理示例来说明该软件包在 HPC 环境中的使用。

2. 动机

大约五年前，普林斯顿大学的地震学小组从事了几个波形反演项目，但应用程序之间的差异和快速的开发速度导致我们的代码变得支离破碎。 这反过来又导致重复工作、易受错误影响以及难以让新成员跟上速度。 需要一个框架来简化我们的软件工作

除了需要灵活的研究工具外，我们还面临着涉及便携性的挑战。 我们的本地集群太小，无法执行我们当时正在工作的所有 3D 反演，并且其他集群的配置与我们自己的不同。 在开始使用新的 HPC 系统时，有时需要花费数周的设置工作才能简单地运行我们的代码。 尽管 SeisFlows 最初是一种研究工具，但用于提高研究灵活性的一些原则也有助于提供跨计算环境的可移植性。

下面，我们将更详细地讨论这些灵活性和可移植性设计目标，然后再进行更全面的封装概述和具有行业基准的实际演示。

2.1. 灵活性

一个可以在不影响可维护性的情况下添加新功能的系统有助于一个包随着时间的推移保持有用，即使在快速发展的研究理念或高流动率的开发团队的情况下也是如此。 涉及模块化设计和面向对象编程的软件工程实践有助于开发灵活的软件（例如，Gamma 等，1995）。 然而，由于遗留代码带来的挑战，这些原则可能难以应用于科学项目。 我们发现，一个可行的回应是开发一个模块化的 Python 框架，其中遗留的可执行文件或原生 Python 代码都可以用于关键组件。 除了提供集成遗留求解器的选项外，Python 还为领域科学家提供了易用性、强大的面向对象功能以及越来越多的科学工具，包括用于线性代数的 numpy (Walt et al., 2011) 和用于地震数据的 obspy 处理（Krischer 等人，2015b）。

2.2. 可移植性

对于 HPC 应用程序，可移植性可能是一个令人生畏的目标，因为集群环境千差万别。 要使 HPC 代码具有可移植性，可能需要能够使用不同的处理器架构、文件系统、内存配置和作业调度程序。

SeisFlows 中的方法是提供一个接口层，工作流通过该接口层与系统资源进行交互。 例如，要启动一组正向模拟，用户通过系统接口的 run 方法调用求解器的正向方法。 通过隔离依赖于环境的属性，系统接口提供了跨不同计算环境的一致命令集。

3. 原理

作为介绍后续章节中使用的术语的一种方式，我们简要回顾一下波形反演（也称为全波形反演或 FWI）的理论。

波形反演是一种数据拟合过程，其中观察数据和合成数据之间的比较用于对对象的内部结构进行成像。 通过观察和合成之间拟合的定量测量，反演问题可以明确地转换为失配函数最小化问题。 misfit 函数的一个常见选择是

$$\chi (m) = \frac{1}{2} \displaystyle \sum_{i=1}^N \displaystyle \int|F[s_i(m;t)] - F[d_i(t)]| ^2 \mathrm{d}{t}, \tag{1}$$

其中 $s$ 是合成数据，$d$ 是观察到的数据，$F$ 是数据处理算子，总和由各种源和接收器 $(i = 1,…,N)$ 我们通常将 $F$ 称为“预处理”运算符，以将其与等式 (2) 中的 $R$ 区分开来。

给定感兴趣对象的材料属性模型 $m$ ，波动方程求解器的任务是通过数值模拟对象内的波传播来生成合成数据 $s$ 参见 Tromp，2015 年的介绍）。

观察到的和合成的数据通常需要某种信号处理来进行有意义的比较。 通常，算子 $F$ 只是一个带通滤波器，但有时需要额外的步骤，例如归一化、静音、去噪或重定数据（有关综合概述，请参阅 Yilmaz，1987 年）。

为了迭代更新模型并改善观察数据和合成数据之间的拟合，我们将非线性优化算法应用于失配函数 $\chi (m)$ 。 每次迭代，都会从之前的模型通过以下方式生成一个新模型

$$
m_ {k + 1} = m_k + \alpha_k H_k^{-1} R[-g_k], \tag{2}
$$
其中 $g$ 是 $\chi(m)$ 的一阶导数的梯度向量，$H$ 是 $\chi(m)$ 的二阶导数的 Hessian 矩阵，或者是依赖于所选非线性优化算法的 Hessian 近似值（参见 Nocedal 和 Wright，2006 年 更详尽的描述）。 最后，$\alpha$ 是由线搜索确定的步长，$R$ 是正则化、平滑或图像处理算子，我们通常将其称为“后处理”算子。

4. 包组织

如图 1 所示，SeisFlows 由六个组件组成。 虽然这种结构是硬连线的，但各个组件本身是高度可定制的

4.1. 求解器

求解器是一种模拟波传播的工具是反演的主要计算引擎。 在模型更新过程中，求解器首先用于生成方程(1)式中的合成数据 $s$ 。 然后通过反投影数据残差（例如，Tromp et al., 2005; Fichtner et al., 2006），求解器用于计算方程(2)中的 $g$ 。

异构介质中的声波或弹性波模拟的选择包括有限差分、有限元和谱元数值方案。 最适合特定反演的求解器类型取决于数据是否涉及体波、表面波或两者兼而有之； 地形或水深变化是小还是大； 以及材料特性在目标结构内是连续变化还是不连续变化。

虽然原则上 SeisFlows 可以与上述所有类型的求解器接口，但迄今为止，主存储库中仅包含谱元求解器接口。 提供了 SPECFEM2D、SPECFEM3D 和 SPECFEM3D_GLOBE 包的接口类，支持各种应用程序。 所有求解器接口都继承自一个通用基类，以避免代码重复。 一些研究人员将 SeisFlows 与有限差分求解器（Gian Matharu，个人通信）结合使用，但此类功能在主包中尚不可用。

4.2. 预处理

在我们的术语中，预处理是指对地震轨迹进行的操作，包括方程（1）中的 $\chi(m)$ 和 $F$ 。 该名称反映了这样一个事实，即此类操作通常在数据残差反投影之前执行。

默认预处理类处理时域数据，并依赖 obspy 来读写数据和信号处理。 如果 obspy 不支持文件格式，用户可以通过简单的插件系统贡献自己的读写实用程序。

默认类中包含各种数据处理选项。 提供高通、低通和带通滤波选项； traceby-trace 或 record 逐个记录部分归一化； 和静音早到或晚到

最后，提供了用于生成伴随迹线的机制，以与“优化然后离散化”求解器一起使用。 提供了对应于波形差、走时、包络和瞬时相位目标函数的伴随轨迹发生器。 可以再次通过插件系统添加对其他目标函数的支持。

4.3. 后处理

在我们的术语中，“后处理”是指对模型、敏感度内核或迁移图像的图像处理操作。该名称反映了此类操作通常在数据反投影之后执行的事实。

各种各样的操作都属于这一类。包括平滑、基础投影、Tikhonov 和总变化正则化在内的技术都涉及灵敏度内核的操作，并且已经使用 SeisFlows 实现（参见 Modrak 和 Tromp，2016 年的概述）。对源或接收器周围的现场效应或数值伪影的处理也属于后处理，在石油和天然气勘探环境中执行的空间过滤或锐化也是如此。

由于后处理操作有时直接在模型或内核上执行，如求解器使用的有限差分、有限元或谱元基础所表示的，求解器和后处理组件之间可能存在一些重叠。在 SeisFlows 中，我们使用 SPECFEM2D、SPECFEM3D 和 SPECFEM3D_GLOBE 包中包含的独立实用程序来平滑基于光谱元素的操作，而不是在 Python 中重新实现这些操作。

4.4. 非线性优化

在生成模型更新和减少数据失配方面，推动反演进展的是非线性优化过程。 波形反演的收敛速度取决于所选择的特定非线性优化算法。

在主程序包中，可用的非线性优化算法包括梯度下降、非线性共轭梯度 (NLCG) 和拟牛顿 (QN) 算法，实施中包含括号和回溯线搜索和重启保护措施（参见 Modrak 和 Tromp，2016 年 概述）。 默认数值设置应该适用于广泛的反演，但也提供了调整参数。 为了提高 NLCG 和 QN 提供的收敛之外的性能，可以通过插件系统加载预处理器。 最后，提供了一个集成测试，该测试使用 Rosenbrock 测试用例 (Rosenbrock, 1960) 检查非线性优化机制。

虽然比完全牛顿法或高斯牛顿法使用的内存要少得多，但共轭梯度和准牛顿算法使用的存储空间大约是简单梯度下降法的 2 到 10 倍。 换言之，NLCG 或 QN 所需的内存是存储单个梯度向量所需内存的 2-10 倍。 尽管求解器操作可以在任意数量的节点上并行执行，但使用 SeisFlows 的非线性优化操作目前仅在单个节点上执行。 如果非线性优化存储成本超过节点上可用的内存，则使用 numpy.memmap 将向量存储在虚拟内存中。 并行化非线性优化例程可以提高非常大的反演的性能，但原因在第 2 节中描述。 4.5，这不是当前的优先事项。

INVERSION WORKFLOW
SOLVER(FORWARD)
PREPROCESS
SOLVER(ADJOINT)
POSTPROCESS
NONLINEAR OPTIMIZATION

MIGRATION WORKFLOW
SOLVER(FORWARD)
PREPROCESS
SOLVER(ADJOINT)
POSTPROCESS

图 2. SeisFlows 提供的默认反演和迁移工作流程。 因为每个组成部分都是高度可定制的，所以上面只显示了一个非常示意性的表示。

4.5. 系统

作为求解器和其他工作流组件与系统资源交互的软件层，“系统”组件提供了从一个环境到另一个环境的互操作性。主包中提供了一系列接口，其中一些用于工作站，另一些用于 LSF、SLURM 和 PBS 集群。在所有情况下，系统组件的想法是提供默认配置，以便在用户环境需要特殊处理的情况下进行定制。此类型自定义的示例可在 github.com/rmodrak/seisflows-hpc 中找到。

以 system.run 函数为例，该函数用于执行令人尴尬的并行任务。在具有单个可用 CPU 内核的计算机上运行时，用户应选择串行系统接口。导致 seisflows.system.serial.run 在运行时加载，它在循环中一次执行一个任务。或者，当在 SLURM 集群上运行小反转时，用户应该选择 slurm_sm 系统接口。这会导致 seisflows.system.slurm_sm.run 被加载，它使用在反转开始时分配的资源同时执行任务。对于较大的反转，用户应该选择以 lg 后缀结尾的系统接口，这会导致资源按需分配，而不是一开始就一次性分配。

特定领域的例子: 全球地震反演
准备球形网格
从数据中删除仪器响应
过滤数据

求解器(波前)

预处理
过滤合成物
选择测量窗口
应用类别权重
比较数据和合成
生成用于“优化-然后-离散化”方法的伴随源

求解器(伴随)

后处理
删除源/接收器工件
合计来源贡献
将 edge_preserving 平滑应用于渐变

非线性优化
初始化拟牛顿求逆 Hessian 矩阵
生成准牛顿搜索方向
执行回溯线搜索

图 3. 可以通过重载默认类在 SeisFlows 中实现的特定领域应用程序的一个示例。 特殊的梯度平滑操作可用于解决全球地震学中地震台站和地震的空间分布不均匀问题。 为此，用户可以从后处理类别中的现有选项中进行选择，或者提供他们自己的自定义后处理类。 上框中列出的步骤仅在反转开始时执行一次。 重复下框中列出的步骤，直到获得所需的失配减少水平。

特定领域的例子: 近地表地震反演
准备笛卡尔求解器网格
筛选数据

求解器(波前)

预处理
过滤合成物
比较数据和合成
生成用于“优化然后离散化”方法的伴随源

求解器(伴随)

后处理
删除源/接收器工件
合计来源贡献
应用 tikhonov 惩罚函数贡献

非线性优化
初始化拟牛顿求逆 Hessian 矩阵
生成准牛顿搜索方向
执行回溯线搜索

图 4. 可以通过重载默认类在 SeisFlows 中实现的域特定应用程序的另一个示例。 针对近地表数据开发了许多专门的数据预处理操作，这些数据缺乏表现良好的地壳、地幔和地核相的全地球数据。 为此，用户可以从预处理类别中的现有选项中进行选择，也可以提供自己的自定义预处理类。

4.6. 工作流

最后，将一切联系在一起的是“工作流程”。提供了默认的反演和偏移工作流程，可以直接使用，也可以作为可以在其上实施专门策略的基类。

在实践中，工作流的执行相当于单步执行包含在 workflow.main 中的代码。用户可以自由修改默认的反演和迁移工作流程，这些工作流程旨在轻松定制。例如，通过包含可用于在模型更新之前或之后进行任何必要设置或清理工作的初始化和完成方法。

图 2 描述了默认的反演和偏移工作流程。因为每个组成部分都是高度可定制的，所以只提供了一个非常示意性的表示。为了了解通过 SeisFlows 可能实现的特定领域应用的类型，图 3 和图 4 描绘了更专业的全地球和近地表反演工作流程。

反映我们小组将 SeisFlows 用作研究工具的情况，通过重载反演和偏移类实现了许多自定义策略。例如，随机反演、双差波形反演和全波形环境噪声反演都可以以这种方式实现。最后，除了执行反演和偏移之外，用户还可以通过以任何所需顺序调用求解器、预处理、后处理和非线性优化组件来编写全新的工作流程。速度分析或不确定性量化任务可以以这种方式实现，举一些例子。

4.7. 杂项

在这里，我们讨论了不属于上述类别的包组织的其他方面。

4.7.1. 测试

Seisflows/tests 包含一组可直接调用的测试脚本。 run_- test_system 使用简单的“hello”消息测试系统接口，run_test_optimize 解决了一个廉价的非线性优化问题，run_test_tools 执行 seisflows/tools 单元测试。

4.7.2. 源代码仓库

代码使用 GitHub (github.com/rmodrak/seisflows) 托管和协作开发。 在拉取请求之后，上述测试使用 Travis 持续集成 (travis-ci.org) 自动运行。 可在 seisflows.readthedocs.org 上查看文档。

4.7.3. 插件

seisflows/plugins 中的每个文件或子目录都包含一组函数，它们共享完全相同的输入和输出语法。 可以通过使用已建立的语法创建新函数来添加新的“插件”。 例如，可以通过模仿现有 misfit 函数的输入和输出，将新的数据 misfit 函数添加到 seisflows/plugins/misfit.py。 因为语法非常简单，所以添加新插件通常需要很少的努力。

4.7.4. 参数和路径

提交工作流时，用户必须指定参数和路径输入文件。 检查值以便可以在不损失队列时间或运行时间的情况下检测到错误。 参数和路径存储在一个字典中，可以在 Python 代码的任何地方访问。

5. 示例：基于 SEAM 山麓模型的陆上油气勘探问题

在过去十年中，陆上地震勘探显着增加，但即使对于最先进的求解器来说，在复杂环境中模拟波传播仍然很困难（Oristaglio，2012 年）。 SEAM II 山麓模型具有山地逆冲带典型的极端地形和强烈的横向和垂直速度变化，是行业联盟为解决陆上建模挑战而创建的三个基准之一（Oristaglio，2013 年）。

下面我们介绍使用山麓数据获得的反演结果。 其中一些结果先前发表在鲍里索夫等人。 （2018）； 在这里，我们关注之前没有讨论过的计算成本、容错性和软件灵活性问题。

5.1. 方法

为了降低实验的总体成本，我们选择了原始模型的 $25 km^2$部分。 所选体积，如图 5 所示，在 x、y 和 z 方向上分别为 $7 \times 3.5 \times 3 km$。 在 1-15 Hz 范围内工作，需要大约 $10^8$ 个数值网格点来进行精确的波模拟。 图 4 显示了模型和数值网格。 高达 0.9 公里的地形变化会导致显着的散射、聚焦和散焦效应。 目标模型包含 P 波和 S 波速度和密度的强烈变化，在某些地方 S 波速度的不连续性超过 1 km/s。 山麓模型的极端地形已被证明会导致有限差分方案出现问题，而诸如 SPECFEM3D 之类的有限元和谱元代码通常表现得更好（Oristaglio，2012 年）。

图5. 近地表地震反演实例。 (a) 目标 S 波速度模型。 (b) 数值网格。 增加了反演的挑战和计算成本，目标模型包含大的地形变化和复杂的折叠和断层结构。 非结构化光谱元素数值网格包含大约 $10^8$ 个积分点。

图6. 近地表地震反演结果 垂直切片（左）和水平切片（右）。 (a,b) 真实模型； (c,d) 初始模型； (e,f) 包络差反演的结果； （ g ， h ）包络差反演的结果，然后是波形差反演。

图7. SEAM II 期山麓模型拍摄记录，位移垂直分量。 左：观察到的痕迹，中：初始合成，右：包络 FWI 后的合成。

在陆上地震勘探中，地表记录的轨迹以高振幅色散地表波为主。与从数据中删除此类信息不同，面波可以用作对浅层结构的附加约束。不幸的是，几个主要障碍，包括面波数值模拟的困难和近地表复杂性的问题影响，使这种方法变得困难。

为了利用面波，反演过程需要灵活性。大多数反演包采用波形差异拟合测量，但这些都不太适合近地表复杂性中的面波反演，加剧了周期跳跃。为了减少循环跳跃，Bozdag 等人提出了包络失配函数。 （2011）。与许多新技术一样，包络目标函数需要在大规模生产之前进行微调。 SeisFlows 为此类测试提供了一个有用的框架。

对于反演，我们使用了 72 个源和 2502 个接收器，它们规则地分布在地表上。接收器之间的距离在 x 和 y 方向上分别为 50 米和 200 分钟，而两个方向上的发射器之间的距离为 600 米。每个镜头对应于在垂直方向上施加的力，其中 Ricker 小波作为源时间函数。我们使用了相对较高的频段（主频为 6 Hz）；我们没有像 Bunks 等人那样使用多尺度方法。 （1995 年）。此外，我们仅对所有接收器使用垂直分量。搜索方向是使用 L-BFGS 算法计算的，该算法在 SeisFlows 中带有提供稳定性的数值保障（例如，Dennis 和 Schnabel，1996）。

尽管裁剪了模型，但反演仍然超出了我们本地分配的容量，因此我们在 Total S.A. 提供的 2000 核 LSF 机器上运行它。只需对 SeisFlows 中的默认 LSF 接口进行少量修改，更改涉及资源队列和 MPI 库路径的规范。网格的内存需求超过了单个节点的容量，因此我们为每个波场模拟使用了多个节点。每个源 6 个节点和每个节点 16 个 CPU 内核，每个源需要 96 个内核。使用 seisflows/system/lsf_lg.py，模拟在作业数组中排队，在资源可用时执行。

5.2. 结果

首先，我们使用包络失配函数来同时反转体波和面波。经过 30 次迭代后，S 波速度模型在浅层地下的精度显着提高（图 6e 和 f）。然后我们将此结果用作波形差异反演的输入，再次使用体波和面波。在进一步的 30 次迭代之后，浅层结构得到了更多的改进（图 6g 和 h）。由于灵敏度集中在自由表面并且振幅比体波大很多倍，因此面波在数据和模型更新方面都占主导地位，因此在面波减弱之前，预计不会恢复深层结构（Borisov 等人，2018 年） .此时使用体波和面波恢复了模型的浅部，可以使用传统的仅体波反演来改善深度恢复。记录位于 x 处的镜头部分 1 公里， 1.75 km如图7所示。由于真实模型的浅部具有很强的非均质性，观测数据以大振幅、色散的瑞利波为主。相比之下，初始合成中的面波的分散性要小得多，因为初始模型相对平滑。很明显，在使用基于包络的反演反演的模型上生成的合成炮集与观察到的更接近。

6. 结论

通过 SEAM II 山麓测试用例，我们说明了 SeisFlows 在一个难题上的使用。 这个例子中的结果部分取决于反演软件的灵活设计，这使得非标准目标函数的实验成为可能。 成功还取决于 HPC 环境之间的可移植性，这使我们能够在超出本地集群分配后切换到 Total S.A. 提供的超级计算机。 除了上述油气勘探基准之外，SeisFlows 的其他用途包括检测非法隧道（Smith 等人，2017 年）； 桥面分层无损检测（Nguyen 和 Modrak，2018 年）； 超声波医学成像（Bachmann，2016）； 利用声学（Modrak 和 Tromp，2015 年）、各向同性弹性模型（Modrak 等人，2016 年）和横向各向异性模型（Rusmanugroho 等人，2017 年）进行油气勘探。 使用 Silwal 和 Tape (2016) 数据的地震层析成像目前正在使用该软件包进行。

致谢

作者感谢 TOTAL S.A. 提供计算资源； SEAM 联盟提供山麓模型； 和 Etienne Bachmann、Wenjie Lei、Gian Matharu、Luan Nguyen、Youyi Ruan、Heru Rusmanugroho、Frederik Simons、James Smith、Yanhua Yuan 和两位匿名审稿人进行了非常有益的讨论。

参考