估计结构对地震运动的响应是一个结构分析中与缓解地震造成的地震风险。许多计算方法结构响应需要了解可以从近地表地震得到的地面研究。在本文中，我们讨论了计算效率波方程层析成像的实现。此方法允许反演初至地震波形进行地震更新可进一步用于评估机械性能的速度模型属性。我们提出了计算效率高的混合用于有限差分（FD）建模的运动学方法初至地震波形。在FD的每个时间步计算仅在移动窄带跟随中执行初至波前。就计算而言，我们得到两个这种方法的优点：计算速度和内存存储计算的首次到达波形时的节省（不是进行计算或存储完整数值所必需的网格）。拟议的方法似乎特别适用于构造广泛用于根据地震数据进行断层速度更新。然后我们应用有效实施初至地震波形的波方程层析成像。
1.简介
地基的力学性质是安全施工的关键信息。特别是，地震分析是高地震风险地区结构分析的重要阶段。这包括估算结构对地震运动的响应。这里的具体主题包括基于性能的抗震设计[1]地震期间隧道和地下结构的动态应力集中、边坡稳定性分析以及输电塔线的设计。许多计算结构响应的方法都需要从近地表地震研究中获得地面力学特性的知识。众所周知，地下地质不规则性可能会严重影响实际场地对地面运动的响应，而地面运动可以模拟复杂的地下结构[2]. 因此，从近地表地震研究中构建详细的地震模型非常重要。
在本文中，我们讨论了用于构建地震速度模型的波方程层析成像的计算效率实现。这种方法需要对地震波场传播（特别是初至波形）进行多次数值计算，这是一个具有计算挑战性的问题，尤其是在3D中。
地震波场建模的数值方法（有限差分、有限元等）在地震数据成像和反演中得到了广泛的应用。它们可以用于计算复杂地下结构中的波场传播，但通常计算成本较高。为此，我们做了大量工作来加速这些方法。
如Boore所述[三]在初至波阵面通过计算网格的某一区域之前，无需计算该区域的波场（在这种情况下，波场均匀为零）。这种称为扩展计算域的方法是由Karrenbach重新提出的[4]. 此外，如果只对首波波形感兴趣，则只需要对首波波前后狭窄区域内的网格点进行数值计算。上述想法由[5]并由Kvasnička和Zahradnřandék进一步开发[6]汉森和雅各布森[7]。
利用有限差分法在规则网格上数值求解航程方程，即可计算出初至波前[8，9]或快速行进法[10]. 注意，与数值求解波动方程相比，计算初至旅行时的成本可以忽略不计。
处理首次到达的类似方法是所谓的SWEET算法，用于使用拉普拉斯域中的阻尼波解来计算首次到达的旅行时间和振幅，如[11]. 该算法可以通过使用核外多波前算法减少每个波长的网格点数量来实现大速度模型[12]. 该算法的优点是可以计算出能量最大的波到达。然而，它与前面提到的方法不同，因为拉普拉斯阻尼会扭曲初至波形。因此，它提供了用于基尔霍夫型成像的有限频率走时和振幅的稳定估计，而不是用于波方程偏移/反演型算法的波形。
计算初到波形的有效方法在构造地震图像的逆时偏移中是有益的[13，14]以及用于建立Luo和Schuster速度模型的波方程层析成像（WT）[15]（传统射线旅行时层析成像的有限频率泛化）。这两种方法都独立地考虑每个炮点集，并意味着计算两个波场：所谓震源场的正向计算和所谓接收场的伴随（时间反转）计算。构造地震图像或速度更新核需要对这两个波场的乘积进行时间积分。因此，需要将前向场的整个历史存储在计算机存储器中，这通常是不可行的。或者，可以在边界处存储前向场的历史，然后与伴随计算同时向后重新计算[16]. 这可以减少用于存储波场的存储器需求，但需要进行额外的计算。对于上述应用，仅在初至波附近的窄带内计算源波场似乎是最佳的，因为它需要的内存少得多。
在本文中，我们首先描述了在一个狭窄的计算频带内计算首波波形的方法，该计算频带跟随首波的运动前沿。然后我们简要回顾了可以利用我们的建模方法的波方程层析成像和逆时偏移方法。最后，我们给出了一些例子，展示了我们在地震成像和速度模型更新中获得的加速。在本文中，我们只考虑二维各向同性模型来说明这一概念。
2.混合正向建模技术
我们进一步发展了在运行窗口中计算波形的思想，并将其应用于层析速度模型的建立，该模型可以用不同的方法完成。最直接的方法是射线层析成像，它只需要对初至波进行旅行时拾取。
关键步骤是按照快速行进法自然产生的时间增长顺序重新排列计算数组[10]. 这种策略允许在开始FD计算之前设置所有计算窗口。因此，我们节省了窗口识别所需的计算机时间，并且在每个FD时间步长时，内存指针都会移到窗口外。此外，使用所提出的策略，可以在每个时间步长将得到的加窗波场存储在RAM中。
2.1、。正向建模的混合方法
让我们列出混合建模方法的主要步骤。为了简单起见，我们考虑了二维波动方程（但很容易将其推广到三维情况和弹性波动方程）。
Vidale提出的窗口加速技术[5]由Kvasnička和Zahradnřandék绘制[6]包括两个步骤：（i）解程函方程得到每个网格节点的旅行时间。（ii）计算第一到达波前周围偏移窗口中波动方程的解，该窗口是使用上一步的旅行时场定义的。
窗口计算技术的数学背景显而易见：影响区域波方程的解由特点.
在每个时间层的第二个波方程求解步骤中，所有网格点都可以更新或保持不变。FD计算仅针对形成初至锋面周围移动区的更新点进行。未更改点中的解决方案取自上一个时间级别。传统方式（Kvasnička和Zahradnřandék[6]、维代尔[5])要确定更新的区域，请检查条件：哪里是时间级别的数量，是FD时间间隔，是第一次到达的旅行时间，以及，是常数。可以使用主波长定义这些常数。
更新条件(1)在波方程FD解算期间，应检查每个时间层的每个网格点，这可能会耗费大量的计算机时间。我们想展示如何避免这种检查，并为基于行程时间递增数组重新排序.
2.2. 数值网格逆风订购（来自Traveltime Solutions）
一如既往，数值解算器提供粘度解与初至旅行时间相对应的eikonal方程。在每个节点一个人有第一次到达的旅行时间。生成的数组可以按递增顺序排序，因此结果将是一个向量（时间递增列表）：
原则上，可以使用任何航程FD解算器或其他计算技术来获得旅行时场用于窗口计算。我们将使用多股快速行进法，由介绍[17]，这是对经典FM的最新且稳健的修改。
这里我们对原作作一个非常简短的描述快速行进法FM每一步的所有网格点集分为三组：接受的分数，窄带、和远点FM算法的主要特点是网格节点中旅行时间计算的特殊顺序，类似于Dijkstra算法。If点是认可的，然后是其中的旅行时间值变得冻结，不再改变。不是的每一点认可的但有一个认可的其邻域中的点应属于窄带; 重新计算这组行程时间。所有其他点是。每次迭代时窄带行程时间的最小值变为认可的.为了提高效率窄带存储在堆结构。
请注意，如果快速行进使用解算器时，给定网格点在时间递增列表中的排名(2)是的号码承兑从FM出来。
让我们考虑一下一对一映射：哪里是一对笛卡尔指数是在旅行时间增加列表中(2);;;.
我们使用映射(三)重新排序所有阵列，用于波方程FD建模。例如，让我们考虑速度模型矩阵; 它可以以向量的形式存储注意映射(三)当需要返回笛卡尔网格时，应记住。例如，可以使用两个整数向量：，。实际上几乎是相同的映射(三)发生在快速行进方法窄带，存储在时间排序的堆结构中。
2.3. 计算和存储孤立波场
对于加窗加速计算，建议的时间增加重排序是自然的。每个时间窗口由其起始编号设置和结束编号使用列表定义的(2). 因此，波方程FD在每个时间层进行计算在循环中完成到无需检查“更新”条件(1).
原则上，基于引入的时间递增重排序，任何时域波方程解算器都可以修改为窗口计算。节点的空间FD模板至少包括一些邻域点（如果使用二阶中心差分）。一旦所有数组重新排序，就应该实现映射(三)多次获得排名用于所有模具节点。更有效的策略是为邻域点的列引入额外的阵列。例如，“左”节点的等级对于节点和军衔可以在数组中记住.
注意，需要为窗口和模具邻域存储一些额外的索引映射数组。另一种方法是，首次波前传播应与计算窗口前面的波方程解同时进行。FM方法提供了这种“动态”旅行时间计算的可能性。然而，应该保持映射(三)（至少向后）。
提出的基于时间递增的阵列重新排序的加窗技术提供了存储波场“历史”的机会。在每个时间级别在当前时间水平上，波场可以逐个记录下来从到在一个大数组中，包含加窗波场的整个“历史”。计算完成后，可以使用存储的寡列极限秩，以时间水平周期读取波场，从而访问波场历史和和使用映射(三)用于获得笛卡尔网格。请注意，可以在相反的时间内生成波场读数。
2.4。进法
考虑各向同性介质中的程函方程：哪里旅行时间在这一点上吗和是地震波速度。我们将在这里简要介绍原作快速行进法在二维情况下。让我们考虑矩形笛卡尔网格：，，、和.原始FM方法[10]基于航程方程的一阶“迎风”FD近似(4):哪里“迎风”近似的本质(5)如图所示1FD模板的网格节点是基于行程时间增加“因果关系”的选择原则；也就是说，“风”的方向就是波前扩展的方向。
FM每一步的所有网格点集分为三组：接受的分数，窄带，和远点FM算法的主要特点是网格节点中旅行时间计算的特殊顺序，类似于Dijkstra算法。If点是接受，那么它的旅行时间值变得冻结，不再改变。不是的每一点认可的但有一个认可的其邻域中的点应属于窄带; 重新计算这一组中的行程。所有其他点是远的。每次迭代时窄带随着旅行时间的最小值认可的.为了提高效率窄带存储在堆结构。
2.5。Virieux交错网格FD方案
本文考虑了二维声波方程：哪里是压力，是密度，以及是波速。
我们使用交错网格[18]解波动方程(7). 为了实现该方案，应将声学方程改写为两个一阶方程：哪里是粒子速度矢量。
考虑的FD方案使用交错网格，因此附加到点笛卡尔网格，附加到“半整数”节点、和压力附加到点此外，速度和压力被附加到不同的时间层。声波方程交错网格FD模板如图所示2FD方程为
3.波方程成像和层析成像的应用
我们展示了所提出的混合计算和存储窗口技术在加速反向时间偏移方面的优势[13，14]和波方程层析成像[15]. 对于这两种方法，我们分别考虑炮点道集，并分别计算两个波场：所谓的震源场是正演计算和所谓伴随场的结果是反向时间计算的结果；即，将时间反转数据作为源签名插入接收器位置。这两个波场的相互作用产生地震图像或速度模型更新。
在每个炮点采集的反向时间偏移中（对应于)我们得到了这两个波场乘积的时间积分图像：
在每个炮点道集的波方程层析成像中（对应于)我们还将灵敏度核构造为两个波场的相互作用：
这里是伴随波场不仅是时间反转数据，而且是时间反转伪残差[15]:哪里是接收器的观测数据，和旅行时间不匹配是由相互相关的观测数据估计的吗和合成数据直接来自波形（最大互相关）：
在总结所有炮点道集的结果时（针对不同的)，我们得到最终的迁移图像：
类似地，我们总结了所有炮点采集灵敏度核（对于不同的)为了获得速度模型更新：
迭代方法可用于速度模型的反演分几个步骤：哪里是处的灵敏度内核第次迭代和正在更新步长。
构建偏移图像的关键步骤(10)或速度模型更新(11)就是计算源场并将其存储起来，以便与伴随场进一步相互关联例如，预计算用于估计用于启动伴随场的伪残差（参见(12)).
因此需要存储前向源字段的整个历史在通常不可行的计算机内存中。或者，可以在边界处存储正向场的历史，然后与伴随计算同时向后重新计算。这可以减少存储波场的内存需求，但需要进行一次额外的波形计算。仅在第一到达附近的窄带中计算源波场对于所提到的应用似乎是最优的，因为它需要相当少的存储器。此外，我们计划仅将建议的方法应用于源场计算伴随接收场时在标准模式下（请注意，它不需要存储在内存中）。
4.示例
4.1. 正向建模
让我们给出一个窗口波场计算的例子。平滑的Marmousi速度模型如图所示三接收机以白色三角形显示，源位置以星星显示。20的Ricker小波 Hz中心频率用作源函数。源1的计算压力数据如图所示4：面板（a）中使用标准FD格式计算的“全”波场，面板（b）中使用提议的混合动力学方法进行“加窗”计算。请注意，第一到达波形看起来是相同的。我们可以看到延迟到达的差异，在“完整”数据上可以看到折射波。对于震源2，计算出的“完整”和“加窗”数据看起来几乎相同，因为在这种情况下没有折射波（因此，我们不显示相应的道集）。
在这个例子中，我们的混合计算方法比标准方法（全波场计算）快20倍。此外，与在每个时间步存储全波场相比，在首次到达后的窄带内存储源波场所需的内存要少17倍。这些加速和内存经济将随着速度模型大小的增加而增加（就主波长而言）。
4.2. 波方程层析成像
让我们考虑由近地表应用驱动的合成示例。在这个例子中，我们使用了梯度模型和高速平方异常，如图所示5（a）。对于这个真实速度，我们使用60 Hz源子波（源和接收器每5个放置一次 m）。对于这些数据，我们使用波方程层析成像进行了速度模型更新（初始模型是具有线性梯度的模型）。使用我们的混合“加窗”方法计算了源波场。速度反演结果(16)如图所示5（b）可以看出，速度异常已恢复。
下一个例子是由地震学中的应用引起的。我们考虑具有两个异常的真梯度速度模型，如图所示6（a）源位置由黑星显示，接收器定期位于表面。合成数据计算为10 Hz里克尔小波。我们使用波方程层析成像进行速度模型更新（初始模型是具有线性梯度的模型）。使用我们的混合“加窗”方法计算了源波场。速度反演结果如图所示6（b）可以看出，速度异常可以恢复。
4.3. 逆向时间迁移
这里我们考虑反向时间迁移示例。计算了Marmousi模型的合成数据，如图所示7（b）; 50次射击，30次 每60个频率放置一个频率为Hz的Ricker源 m。偏移速度模型如图所示7（a）.
我们使用“窗口”混合方法来建模源字段; “窗口”建模比“完整”波形建模快33倍（在Intel Core i7-4770 3.4上分别为1.5秒和50秒 GHz CPU）。它还需要20倍的内存来存储首次到达的波场，即0.8 Gb与17 千兆字节。50个炮点道集的反向时间偏移图像(15)如图所示8：使用“完整”波形计算（a）和使用我们的“加窗”计算（b）。这两个迁移结果看起来完全相同。注意接收器波场在这两种情况下都是以“完整”波形模式计算的。
5.讨论
示例表明，我们的“窗口化”方法允许将任何2D模型的“窗口”源波场存储在内存中，而在大多数情况下，在内存中拟合“完整”源波电场是不可行的。使用我们的“加窗”方法进行波方程层析成像和逆时偏移的总体速度接近3。对于标准方法，源场不适合内存，因此需要三个标准的“完整”波场计算[16]用于实现层析成像更新或迁移。我们的“窗口化”方法大大减少了源场计算所需的时间，并且可以将其存储在内存中以供进一步使用。因此，我们的“加窗”方法的计算成本接近于接收器场的一个“完整”波场计算。
我们想提到在对比模型中应用“窗口化”建模方法的潜在问题[6]. eikonal解算器提供第一次到达，这可能对应于弱头波，而用于成像的高能波可能稍后到达。第一次到达锋面之后的计算窗口可能会错过这些波。然而，我们的方法所适用的速度模型的类别足够大，可以用于实际应用。
6.结论
我们提出了一种新的混合动力学有限差分（FD）方法，用于计算和存储初至地震波形。在这种“窗口化”的计算方法中，在每个时间步长，我们只在第一个到达波前之后的狭窄区域中执行FD计算。该波前是使用快速行进的eikonal解算器预计算的。我们提出了一种新的有效的计算和波场存储策略，该策略基于时间递增的重新排序。
我们展示了在波方程层析成像和偏移中成功使用“窗口”建模方法的几个示例。与标准的“全”波场计算和存储相比，该方法的计算速度提高了20–30倍，存储源波场的内存减少了17–20倍。实现波方程层析成像或逆时偏移的总速度提高了3倍。
在本文中，我们只考虑二维各向同性模型来说明这一概念。扩展到3D模型很简单。各向异性模型方法的推广需要有效的各向异性eikonal解算器算法。这是我们未来研究的主题。
利益冲突
作者声明，本论文的出版不存在利益冲突。
致谢
研究得到了俄罗斯基础研究基金会（14-05-00862号拨款）的部分支持。作者非常感谢这位不知名的评论员的评论，这些评论帮助改进了论文。
