A frequency domain elastic wave full waveform inversion method based on nearly analytic discrete operator
-
摘要:
全波形反演是一种利用地震波传播的动力学特征来获取地下介质物性参数的反演方法，可为揭示地下精细结构提供重要依据。本文以弹性波方程作为数学模型来模拟地震波传播规律并进行相应的反演方法研究。为提高计算效率与反演结果的准确性，可将近似解析离散化（NAD）算子用于频率域弹性波方程的正演模拟。本文在频率域NAD离散的基础上推导阻抗矩阵的稀疏分块结构与反演目标函数对模型参数的梯度计算公式，由此建立基于NAD算子的频率域弹性波全波形反演方法。为验证该方法的有效性，文中通过数值实验对多种典型介质模型进行反演计算，均得到了理想的反演结果。
Abstract:Since ancient times, earthquakes have been an unavoidable natural phenomenon. As one of the regions with frequent seismic activities, China urgently needs advanced underground detection technology to characterize the fine structure of the lithosphere. This will provide a better understanding of the generation and development laws of earthquakes and minimize their impact on human society. Full waveform inversion （FWI） is an inversion method that utilizes the dynamic characteristics of seismic wave propagation to obtain subsurface physical parameters, providing important insights into the fine structure of the subsurface. This method further enriches the theoretical framework of conventional velocity modeling methods and has gradually become a research hotspot due to its high accuracy, multi-parameter, and multi-dimensional modeling advantages. As a high-resolution seismic inversion method, FWI minimizes the objective function that measures the discrepancy between the synthetic waveform and seismic observation data through wavefield simulation and optimization iteration, obtaining optimal subsurface medium parameters and constructing imaging. In the early stages of research, FWI was primarily performed in the time domain. Subsequently, researchers extended the inversion theory to the frequency domain. In comparison, frequency-domain inversion methods are conducive to parallel computing, reducing computational burden, and easy characterization of attenuation effects. By selecting a small number of frequencies that sufficiently cover the wavenumber, the fine structural details of the inversion model can be characterized without compromising inversion accuracy, thereby further reducing the computational cost of frequency-domain FWI. Additionally, frequency-domain wavefield simulation exhibits higher efficiency, especially in the case of multiple sources and increased attenuation factors. Therefore, conducting FWI in the frequency domain holds significant importance. With the continuous development of seismic research and rapid improvement in computer hardware, traditional acoustic wave equations have difficulty meeting the requirements of fine subsurface detection. FWI has gradually evolved towards the direction of elastic parameter inversion based on isotropic, anisotropic, viscoelastic wave equations, among others. Compared to acoustic waves, elastic waves include the propagation laws of P-waves and S-waves, providing richer wavefield information. Furthermore, in the near-surface region where the propagation time is short, it is challenging to separate body waves and surface waves, and attenuating or filtering surface waves is not easy. Both types of waves need to be considered for imaging. The elastic wave equation accurately simulates seismic wave propagation in subsurface media. Therefore, this study adopts the elastic wave equation as a mathematical model to simulate the propagation laws of seismic waves and conducts corresponding research on inversion methods. In FWI, the main computational burden lies in the forward modeling process, and the computational efficiency and accuracy of the inversion heavily depend on the quality of the forward simulation method. Considering the complex structure of the elastic wave equation, traditional numerical solving methods require more computational costs. Therefore, employing efficient forward simulation methods is crucial for improving computational efficiency and the accuracy of inversion results. The nearly analytie discrete （NAD） operator has emerged as an important class of differencing operators in recent years. Its construction principle involves approximating high-order partial derivatives using the cross-combination of nodal displacements and their gradient values. The NAD operator offers high accuracy and significantly enhances computational efficiency. Therefore, in order to improve computational efficiency and the accuracy of inversion results, this study introduces the NAD operator for frequency-domain elastic wave equation forward simulation. This study establishes a frequency-domain elastic wave FWI framework based on the NAD operator. It derives the gradient calculation formula for the inversion objective function using the sparse block structure of the impedance matrix. The idea of separating compressional and shear waves is implemented throughout the entire solution process, greatly reducing the scale of gradient calculation and suppressing interference between them, further improving inversion computational efficiency. Finally, numerical experiments using three classic media models are conducted, obtaining ideal inversion results. The effectiveness and correctness of the proposed elastic wave inversion method are validated through the trend curves of inversion errors at different frequencies. In summary, the frequency-domain elastic wave FWI method based on the NAD operator derived in this study exhibits significant advantages in the inversion of subsurface physical parameters. Through numerical experiments, we have demonstrated the effectiveness and accuracy of this method in capturing fine subsurface structures. We believe that this method has broad application prospects in geological exploration and resource development and provides valuable references for researchers in related fields.
-
引言
随着我国数字地震观测技术的快速发展，台站密度逐渐增大，地震监测能力和水平得到大幅度提高，地震台网能够捕获更多事件信号。既可以记录到震级较小的天然地震事件，又可以记录到大量非天然地震事件（赵永等，1995），如会导致地震预警误触发的人工爆破等。及时准确地从天然地震目录中剔除人工爆破事件以保证目录的完备性和准确性对于区域强震预测及地震危险性评估至关重要（黄汉明等，2010；隗永刚等，2019）；快速有效且自动区分天然地震和人工爆破，对防震救灾工作中的抢险应急响应具有重要意义，已成为目前地震监测工作中面临的重要问题之一，也是近年来震灾救援部门和地震学者高度重视的问题。
随着AlphaGo和ChatGPT风靡全球，让人类重新思考人工智能的更多发展可能性，对于人类体力劳动和非创造性脑力劳动，利用人工智能的手段替代能很大程度的提高效率、节省人力成本并减少错误。例如在地震学上，研究学者应用人工智能手段进行地震波降噪、地震震源机制估算、微震识别及定位等（Gao，Zhang，2019；Kuang et al，2021；Zhang et al，2022）。近年来，部分学者尝试利用人工智能技术进行地震波形分类。例如，应用神经网络模型对地震和噪声波形进行分类，准确率普遍达到95%以上（赵明等，2019；Chen et al，2019；Men-Andrin et al，2019；Zhang et al，2020a，b）。其中，Men-Andrin等（2019）应用了五种神经网络模型对实时地震波形和噪声进行判别，卷积神经网络（convolutional neural network，缩写为CNN）的表现最为突出，拥有最高的精确率（precision）和召回率（recall）。不少学者也应用卷积神经网络模型进行天然地震与人工爆破的识别，例如：陈润航等（2018）从震源波形中提取梅尔频率倒谱系数图后用卷积神经网络进行天然地震与爆破事件分类识别，准确率（accuracy）达95%以上；Linville等（2019）利用卷积神经网络模型对美国犹他州的天然地震信号和采石场爆破信号进行区分，区分结果可达99%；周少辉等（2021）利用深度学习技术中的卷积神经网络模型对山东地区的天然地震与非天然地震进行学习训练，训练验证识别准确率达到98%以上；郑周等（2023）利用卷积神经网络对福建及邻省的地震、噪声、爆破及异常波形进行区分，识别率分别达到97.9%，99%，99.2%和99.3%。上述工作都证明深度学习在地震学上存在一定的发展潜力，卷积神经网络在波形区分上的表现更为突出。
目前广东省地震台网一周平均记录到爆破事件72次，且爆破识别工作完全通过人工完成，不仅耗时而且还高度依赖分析人员的专业水平，难免会出现错误，如广东丰顺地区的地震波形与爆破类似，极易判断错误。再者，随着预警系统的正式运行，全省预警速度达到秒级，需要我们对台站收集的各类信号进行快速并自动分类，以避免发生误触发。由于每个地区、省份的地质构造不同，其天然地震波形甚至包括爆破的其它非天然地震事件波形也大相径庭，即使是具有高泛化能力的卷积神经网络也不太可能将同一套训练好的模型用于区域或更大的范围。
为解决以上问题，本文拟搭建一个适用于广东地区的卷积神经网络模型，利用广东省地震台网的天然地震与人工爆破数据训练该模型，构建一个适用于广东地区的人工爆破智能识别器，以期简化数据预处理流程，缩短事件判定时间，高效、准确、稳定判别广东地区天然地震与人工爆破波形。
1. 数据资料
由于广东省人工爆破工程作业的日益增多，人工爆破成为广东省地震台网非天然地震事件记录最多的事件类型。2013—2023年广东省地震台网共记录到ML＞1.8的人工爆破事件472次， 其中发生在省内的事件达312次（图1红色圆圈）。而广东测震台网观测台站共112个，其中包含41个省外共享台，每个台站记录的波形均包括垂直、东西、南北三分向，台站观测数据采样率均为100 Hz。本文从广东测震台网选取2013—2023年发生在广东省内的人工爆破事件，并为平衡天然地震和人工爆破事件的数量，选取2017年6月至2019年1月期间发生在广东省陆地内ML＞1.4的526次天然地震事件（图1蓝色圆圈）。这些地震事件均由广东省地震局台网中心人员进行过专业的地震编目分析并入库。由于人工爆破面波能量衰减很快，为了达到更好的识别效果，选取震中距30 km内的事件波形。
![]()
图 1 所选事件空间分布图及广东省内台站分布图Figure 1. Spatial distribution of events uesd in this paper and seismic stations in Guangdong研究样本根据事件P波到时信息截取：截取P波到时前5 s的噪声记录、总长为60 s的波形，并能完整包含P波、S波，在事件资料充足的情况下，为满足快速识别的要求、尽可能减少识别时间以及达到更好的识别效果，最终选取单台单向（垂直向）波形为训练数据。去除标定、方波、突跳、干扰和仪器故障等异常信息后，共得到
1840 份有效波形，包括天然地震事件1 446份、人工爆破394份。选取其中1 000份天然地震事件波形和300份爆破波形用于训练学习，组成训练集；剩余用于测试训练后的模型，组成测试集。图2给出了天然地震（图2a）与人工爆破（图2b）的波形样本图。![]()
图 2 天然地震（a）和人工爆破（b）的有效波形示例Figure 2. Examples of effective waveforms for earthquake （a） and blasting （b）利用神经网络较强的分类能力，将本研究内容变为天然地震与人工爆破的二分类问题，并建立样本标签。天然地震作为第一类，输出标签结果为0；人工爆破为第二类，输出标签结果为1。
2. 卷积神经网络模型
2.1 卷积神经网络的原理
Hinton等（2006）最早提出深度学习的概念，意指基于样本数据通过一定的训练方法得到包含多个层级的深度网络结构的机器学习过程（Bengio et al，2013）。深度学习所得到的深度网络结构类似于神经网络，故称之为深度神经网络。深度神经网络分为三类：前馈深度网络、反馈深度网络以及双向深度网络。典型的前馈深度网络包括多层感知机、卷积神经网络等。
卷积神经网络可以用来处理具有类似网格结构的数据（LeCun et al，2015），不仅容错性好、自适应性与自学能力强，而且还可将图像作为数据输入，具有自动识别图像特征、参数估计数量较少、可以采用权值共享网络结构来降低模型复杂度等优势，提高了图像识别效率和准确率。Krizhevsky等（2012）首次将卷积神经网络应用于ImageNet大规模视觉识别挑战赛（imageNet large scale visual recognition challenge，缩写为ILSVRC），并取得图像分类和目标定位任务第一名。
CNN以原始数据作为输入，通过卷积、池化和非线性激活函数映射等一系列计算，将原始数据逐层抽象为自身任务所需的最终特征表示，最后以特征到任务目标的映射作为结束。卷积神经网络的基本结构包括输入层（input）、卷积层（convolutional layers）、池化层（pooling layers）、全连接层（fully-connected layers）和输出层（output）（图3）。卷积层由多个滤波器组成，进行卷积计算，旨在学习输入的特征表示，并通过激活函数传递给下一层。激活函数又称非线性映射函数，非线性被引入后实现了多层网络对非线性特征的检测，典型的激活函数有sigmoid，tanh和ReLU等（Han，Moraga，1995；Glorot，Bengio，2010；Nair，Hinton，2010）。池化层位于两个卷积层之间，其作用是将卷积层提取的事件特征信号进行降维处理，防止出现过拟合现象。典型的池化操作有最大池化、平均池化和随机池化。全连接层是当前层每个神经元与上一层所有神经元的连接，目的在于产生全局语义信息。卷积层的输出特征公式为：
![]()
图 3 本文采用的卷积神经网络结构Figure 3. The architecture of the proposed convolutional neural network$$ x_j^l ＝ f\left(\sum\limits_{i \in {M_j}} {{x^{l － 1}}k_{ij}^l ＋ b_j^l} \right) \text{，} $$ (1) 式中，$ {M_j} $为上层特征子集，$ l $为层数，$ x_j^l $为第$ l $层的第$ j $个特征图，$ b $表示偏置，$ k $表示卷积核，$ f $表示激活函数。
AlexNet卷积神经网络模型较早被提出且至今仍被广泛应用，该网络模型为2012年Krizhevsky赢得ILSVRC比赛冠军所用的模型，其主要由5个卷积层和3个全连接层组成，共包含约65万个神经元以及6 000万个可训练参数（Krizhevsky et al，2012）。为了适应小批量的训练集和验证集数据，更好地适应广东地区地震事件特征分类，本文将AlexNet模型的最后三层替换为全连接层、softmax层和分类输出层，并在全连接层中设置分类数目及增大其学习率因子以达到加快最终层学习速度的目的。Softmax函数将分类输出中的概率进行归一化处理，计算式为：
$${q_i} （ x ） ＝ \frac{{{{{\mathrm{e}}} ^{{{\textit{z}}_i} （ x ） }}}}{{{{{\mathrm{e}}} ^{{{\textit{z}}_1} （ x ） }} ＋ {{{\mathrm{e}}} ^{{{\textit{z}}_2} （ x ） }}}}\text{，} $$ (2) 式中，$ i ＝ 0\text{，} 1 $分别表示天然地震事件和人工爆破事件，$ {{\textit{z}}_i} （ x ） $是最后一层全连接层输出的未归一化的$ i $类概率。两个卷积层间采用的激活函数为ReLU （Nair，Hinton， 2010），其将特征映射中的所有负像素值替换为零的非线性操作极大地减少了计算量，加快了收敛速度；池化层采用最大池化函数（Zhou，Chellappa，1988），池化层各输出特征图计算方法为：
$$ x_j^l ＝ f ［ \beta _j^iD （ x{}_j^{l - 1} ） ＋ b_j^l ］ \text{，} $$ (3) 式中，$D （ ） $表示池化函数，每一个输出特征图对应一个权重系数$ \beta $和偏置$ b $。图3为本文所采用的卷积神经网络结构示意图。神经网在双GPU上运算，输入图片大小为227×227×3；各卷积层采用的卷积核大小依次为11×11，5×5，3×3，3×3，3×3；各池化层采用的最大池化操作都是3×3，即以3×3的窗口大小进行缩放，获取不同粗糙程度的特征。
2.2 代价函数
模型训练过程中，为了让实际网络输出更准确，需要调试各层的权值，故使用交叉熵作为代价函数计算预测值与真实值之间的误差并反向传播给各层，并更新各层权值，其计算方法为：
$$ H （ p\text{，} q ） ＝ - {\sum\limits_x} {p （ x ） \ln q （ x ） } \text{，} $$ (4) 式中：$ p $为真实概率值；$ q $为预测值，由CNN最后一层输出值经过Softmax函数计算后所得。
2.3 搭建卷积神经网络模型
针对波形分类问题，采用单台单通道60 s波形作为输入，输出判别结果为0 （天然地震）或1 （人工爆破）。卷积神经网络利用卷积操作提取波形特征，经过一系列计算和处理，最后采用类似投票的形式输出对应波形类别的概率。基于AlexNet卷积神经网络模型搭建的快速高效识别天然地震与人工爆破的识别器主要由两个模块组成：训练模块和测试模块，图4给出了两个模块的简化图。
![]()
图 4 基于AlexNet卷积神经网络模型搭建的人工爆破识别器流程图Figure 4. AlexNet-based blasting classifier flow chart训练模块的核心步骤包括：
1） 在模型构建之前，挑选信噪比高、记录清晰、波形震相明显且两类事件特征较突出的数据作为训练学习的素材，这将有利于模型快速抓取事件特征、快速生成适合区分广东省内地震波型类型的参数集。
2） 为训练数据集添加相应的事件类型标签。
3） 将用于训练的各类事件波形作为输入，利用AlexNet卷积神经网络对事件进行识别训练和验证。前人的研究表明概率阈值的设定对事件判别的准确率影响不大（Li et al，2018；郑周等，2023），即CNN对事件判定较为准确：对一个图片的判断概率值基本接近1或0，中间值很少。因此，对于本文研究的二分类，将概率阈值设为0.5。单次训练周期中根据迭代的次数随机选择计算验证集代价函数以及准确率，当验证集代价函数曲线趋于直线且识别准确率几乎不变时停止训练。随后根据验证集呈现的结果进行效能评估和模型优化，包括改变初始学习率、调整迭代次数、改变梯度下降速度等，使该模型更好地适应广东地区非天然地震的识别任务，进而让识别结果达到最优。
将用于测试的事件波形随机打乱后输入经参数优化后的卷积神经网络模型中，根据模型输出的事件概率大小，判定波形为天然地震事件或者人工爆破，构建一个自动识别器。
3. 结果分析
3.1 评价指标
虽然准确率A可以判断总的正确率，但在样本不均衡的情况下，并不能作为衡量结果的最优指标；在文中天然地震事件与人工爆破事件样本不均衡的情况下，引入精确率P、召回率R和F1分数（F1-score）对模型结果进行全面评价。其中精确率表示模型正确预测为正的样本的数量；召回率表示正例样本被预测出的数量；F1分数则是对精确率和召回率的一种调和平均。计算式为：
$$ A ＝ \frac{{{\mathrm{TP}} ＋ {\mathrm{TN}}}}{{{\mathrm{TP}} ＋ {\mathrm{TN}} ＋ {\mathrm{FP}} ＋ {\mathrm{FN}}}} \text{，} $$ (5) $$ P ＝ \frac{{{\mathrm{TP}}}}{{{\mathrm{TP}} ＋ {\mathrm{FP}}}}\text{，} $$ (6) $$ R ＝ \frac{{{\mathrm{TP}}}}{{{\mathrm{TP}} ＋ {\mathrm{FN}}}}\text{，} $$ (7) $${ F_1 }＝ {\frac{{2PR}}{P＋R}} \text{，} $$ (8) 式中：$ {\mathrm{TP}} $和$ {\mathrm{FP}} $分别代表真阳性（true positive，即天然地震波形正确预测为天然地震波形）和假阳性（false positive，即人工爆破波形错误预测为天然地震波形），而$ {\mathrm{TN}} $和$ {\mathrm{FN}} $代表真阴性（true negative，即人工爆破波形正确预测为人工爆破波形）和假阴性（false negative，即天然地震波形错误预测为人工爆破波形）。
3.2 验证集比例及训练样本数量对准确率的影响
在模型的训练模块中，将输入的训练数据集分为训练集和验证集，两者常见比例为7 ∶ 3或8 ∶ 2。为了得到最好的训练效果，先按不同比例分割训练所用数据，再进行测试，使用训练后验证集的准确率来评估效果。同样，为了测试不同训练样本数量对训练效果的影响，改变数据量比较大的天然地震事件的训练数量，查看验证集准确率的改变。
如图5a所示：当训练集与验证集比例大于3 ∶ 7后，根据式（5）计算的准确率A随比例小幅度增大，人工爆破的准确率在波动但总的趋势为随比例增大而增大，且天然地震事件的准确率高于95%；当比例为8 ∶ 2时，总准确率达到最大97.31%；当比例大于8 ∶ 2时（9 ∶ 1），总准确率呈下降趋势。可见，训练过程提供给模型学习的样本数应多于其用于验证模型参数以及评估模型效果的数量，此模型中训练集与验证集比例为8 ∶ 2时准确率达到最佳。如图5b所示，随着用于训练的天然地震事件数量增加，其准确率增加，当数量大于600后，准确率高于95%。人工爆破事件的准确率虽随地震数量增加而下降，但其准确率始终高于90%，如果有数量足够多的爆破素材且与天然地震数量相当，预计其准确率有望达到95%以上，模型学习效果也会更佳。
![]()
图 5 训练过程中训练集与验证集比例（a）及天然地震训练数量（b）与验证集准确率的对应关系Figure 5. Validation accuracies verus ratios of training and validation （a） and validation accuracies verus training sizes of natural earthquakes （b）3.3 神经网络训练结果
将用于训练学习的
1000 份天然地震波形和300份爆破波形按照8 ∶ 2的比例随机分为训练集和验证集后导入AlexNet卷积神经网络模型，详细训练过程如图6所示，整个训练过程用时7分4秒。由图可知，随着训练次数的增加，训练集的准确率逐渐上升并超越90%，最终稳定在某一数值附近；与此同时，代价函数曲线快速下降并最终稳定在相对较小的数值附近不再变化；需要强调的是，模型未发生过度拟合现象。最终，训练集准确率高达100%，验证集为97.31%。![]()
图 6 AlexNet卷积神经网络模型训练结果（a） 准确率；（b） 代价函数损失Figure 6. Training performance of convolutional neural network of AlexNet（a） Accuracy；（b） Cost function loss3.4 测试结果
对训练好的模型进行测试，将剩余的446份天然地震与94份人工爆破事件的波形随机打乱后导入模型中进行事件类型识别，识别完成花费时间不到2 s，准确识别样本数526个，准确率达到97.41%。用正确分类到不同波形类型的百分比来评估该网络模型的鲁棒性，从表1可以看出：天然地震的精确率达到98%以上，爆破精确率达到90%以上，得出该模型能有效识别广东地区的天然地震与人工爆破波形。
表 1 测试得到的地震与爆破的精确率、召回率和F1分数Table 1. Precision，recall and F1-score for two classification categories事件类型 精确率 召回率 F1分数 天然地震 0.989 0.980 0.984 人工爆破 0.908 0.947 0.927 如混淆矩阵（表2）所示（每一行之和表示该类别的真实样本数量，每一列之和表示被预测为该类别的样本数），天然地震波形出现了9次分类错误，人工爆破出现了5次分类错误，从比例可以看出，该模型对天然地震波形的识别率优于人工爆破。
表 2 测试识别结果的混淆矩阵Table 2. Confusion matrix of waveform classification天然地震预测值 人工爆破预测值 天然地震真实值 437 9 人工爆破真实值 5 89 4. 讨论与结论
为解决广东台网目前费时又考验工作人员专业知识的靠人工的事件分类难题，并且为探索卷积神经网络模型在天然地震与人工爆破识别应用中的泛化性，本文构建了AlexNet卷积神经网络模型，并利用广东省台网编目入库事件震中距30 km内的单台垂直向波形，训练模型并优化参数后进一步测试该模型，尝试搭建人工爆破自动识别器，获得如下认识：
1） 对训练好的模型进行测试，共 540份样本中，正确识别526份，准确率达97.41%，其中天然地震事件和人工爆破事件的精确率、召回率以及F1分数分别均大于0.98和0.90。从实际应用而言，因为广东地区的爆破震级不大，所以保证地震目录完备性和准确性更具重要意义，因此要求天然地震事件的精确率和召回率高于爆破，实验结果说明AlexNet卷积神经网络模型对地震与爆破的识别率高，且拥有基本不漏地震事件的能力，具有较好的实际应用价值；
2） 训练AlexNet卷积神经网络模型所需要的样本不仅前期处理简单且所需训练的数据量少，当样本数多于600时，准确率可高于95%，虽然此处的样本数量600不一定存在普适性，但可以看出训练样本数达到百级后，该卷积神经网络便可得到较高的准确率，颇具高效性；
3） 随着天然地震样本数增加，天然地震识别率上升，人工爆破下降（图5b），这是因为模型接收了更多天然地震的信息和特征，对地震的适应性更高，随着时间推移，爆破素材不断累积，便可持续不断往模型中输入新的训练数据，这将不断提高爆破事件识别的适应性和准确性，进一步提高其识别率。
本文模型训练所用人工爆破事件数量不够充足，该实验结果并不能代表卷积神经网络模型未来的实际应用效果。但实验结果显示：基于AlexNet卷积神经网络模型对广东天然地震与人工爆破的事件波形类型识别准确率达97%以上，天然地震事件的精确率、召回率以及F1分数均大于0.98，能够满足及时准确地从天然地震目录中剔除人工爆破事件保证目录完备性和准确性的要求，有利于区域强震预测及地震危险性评估；另外，随着后续更多爆破事件素材的输入，将进一步提高模型对人工爆破的识别率，基于该模型搭建的爆破识别器将能为广东省地震台网工作节省大量时间和人力成本，并为震后应急工作提供类型结果支撑。未来研究方向将围绕识别器实用化展开，立足于广东省地震台网记录资料，通过不断检验，进一步提高识别器的准确性和稳健性，并将其应用于实际地震预警与日常地震监测中。
-
![]()
图 15 P波速度vP （左）和S波速度vS （右） Marmousi模型第一阶段（a）和第二阶段（b）反演结果图
Figure 15. The first-stage （a） and the second stage （b） inversion results of Marmousi model for P-wave velocity vP （left） and S-wave velocity vS （right）
![]()
图 1 频率域NAD全波形反演流程图
Figure 1. Flow chart of NAD full waveform inversion in frequency-domain
![]()
图 2 P波速度vP （a）和S波速度vS （b）异常块反演真实模型
Figure 2. Real model of anomaly block inversion for P-wave velocity vP （a） and S-wave velocity vS （b）
![]()
图 4 P波速度vP （左）和S波速度vS （右）异常块模型第一阶段（a）及第二阶段（b）反演结果图
Figure 4. The first-stage （a） and the second-stage （b） inversion results of the anomaly block model for P-wave velocity vP （left） and S-wave velocity vS （right）
![]()
图 6 P波速度vP （a）和S波速度vS （b）异常块模型反演结果图
Figure 6. The inversion results of the anomaly block model for P-wave velocity vP （a） and S-wave velocity vS （b）
![]()
图 7 P波速度vP （a）和S波速度vS （b）双层介质真实模型
Figure 7. Real model of double-layer medium for P-wave velocity vP （a） and S-wave velocity vS （b）
![]()
图 9 P波速度vP （左）和S波速度vS （右）双层介质模型第一阶段（a）和第二阶段（b）反演结果图
Figure 9. The first-stage （a） and the second stage （b） inversion results of the two-layer medium model for P-wave velocity vP （left） and S-wave velocity vS （right）
![]()
图 12 P波速度vP （a）和S波速度vS （b）双层介质模型反演结果图
Figure 12. The inversion results of the two-layer medium model for P-wave velocity vP （a） and S-wave velocity vS （b）
![]()
图 13 P波速度vP （左）和S波速度vS （右） Marmousi真实模型（a）和初始模型（b）
Figure 13. Marmousi real model （a） and initial model （b） for P-wave velocity vP （left） and S-wave velocity vS （right）
![]()
图 17 P波速度vP （a）和S波速度vS （b） Marmousi模型反演结果图
Figure 17. The inversion results of Marmousi model for P-wave velocity vP （a） and S-wave velocity vS （b）
表 1 三种介质模型的反演参数设置
Table 1 Inversion parameter settings of three media models
正演网格点数（nx×nz） 空间步长h/km 计算区域/km2 介质密度ρ/（kg·m−3） 异常块模型 101×81 0.025 1.75×2.25 2.63 双层介质模型 101×81 0.025 1.75×2.25 2.63 Marmousi模型 233×78 0.015 5.775×1.9 2.63 
下载: 导出CSV
-
卞爱飞,於文辉,周华伟. 2010. 频率域全波形反演方法研究进展[J]. 地球物理学进展,25(3):982–993. doi: 10.3969/j.issn.1004-2903.2010.03.037 Bian A F,Yu W H,Zhou H W. 2010. Progress in the frequency-domain full waveform inversion method[J]. Progress in Geophysics,25(3):982–993 (in Chinese).
黄少华,任志明,李振春,谷丙洛,李红梅. 2019. 纵横波分离的多震源弹性波全波形反演[J]. 石油地球物理勘探,54(5):1084–1093. Huang S H,Ren Z M,Li Z C,Gu B L,Li H M. 2019. Multi-source elastic full waveform inversion based on P-wave and S-wave separation[J]. Oil Geophysical Prospecting,54(5):1084–1093 (in Chinese).
郎超,仇楚钧,刘少林,申文豪,李小凡,徐锡伟. 2021. 求解弹性波动方程的频率域近似解析离散化波场模拟方法[J]. 地球物理学报,64(8):2838–2857. doi: 10.6038/cjg2021O0332 Lang C,Qiu C J,Liu S L,Shen W H,Li X F,Xu X W. 2021. A nearly discrete analytic method of wave-field simulation for elastic wave equations in the frequency domain[J]. Chinese Journal of Geophysics,64(8):2838–2857 (in Chinese).
郎超,刘少林,杨晓婷,徐锡伟. 2022. 基于改进NAD方法的频率域声波逆时偏移[J]. 地球物理学报,65(3):1071–1085. doi: 10.6038/cjg2022P0501 Lang C,Liu S L,Yang X T,Xu X W. 2022. Frequency-domain acoustic reverse time migration based on improved NAD method[J]. Chinese Journal of Geophysics,65(3):1071–1085 (in Chinese).
刘璐,刘洪,刘红伟. 2013. 优化15点频率-空间域有限差分正演模拟[J]. 地球物理学报,56(2):644–652. Liu L,Liu H,Liu H W. 2013. Optimal 15-point finite difference forward modeling in frequency-space domain[J]. Chinese Journal of Geophysics,56(2):644–652 (in Chinese).
龙锋,闻学泽,徐锡伟. 2006. 华北地区地震活断层的震级-破裂长度、破裂面积的经验关系[J]. 地震地质,28(4):511–535. doi: 10.3969/j.issn.0253-4967.2006.04.001 Long F,Wen X Z,Xu X W. 2006. Empirical relationships between magnitude and rupture length,and rupture area,for seismogenic active faults in North China[J]. Seismology and Geology,28(4):511–535 (in Chinese).
齐诚,赵大鹏,陈颙,陈棋福,王宝善. 2006. 首都圈地区地壳P波和S波三维速度结构及其与大地震的关系[J]. 地球物理学报,49(3):805–815. doi: 10.3321/j.issn:0001-5733.2006.03.024 Qi C,Zhao D P,Chen Y,Chen Q F,Wang B S. 2006. 3-D P and S wave velocity structures and their relationship to strong earthquakes in the Chinese Capital region[J]. Chinese Journal of Geophysics,49(3):805–815 (in Chinese).
任志明. 2016. 声波和弹性波波动方程有限差分正反演方法研究[D]. 北京:中国石油大学(北京):2−3. Ren Z M. 2016. Research on Finite-Difference Modeling and Inversion Methods Based on Acoustic and Elastic Wave Equations[D]. Beijing:China University of Petroleum (Beijing):2−3 (in Chinese).
宋海斌,张关泉. 1998. 层状介质弹性参数反演问题研究综述[J]. 地球物理学进展,13(4):67–78. Song H B,Zhang G Q. 1998. A summary of inversion problem research on elastic parameters of stratified media[J]. Progress in Geophysics,13(4):67–78 (in Chinese).
孙成禹,李振春. 2011. 地震波动力学基础[M]. 北京:石油工业出版社:1−2. Sun C Y,Li Z C. 2011. Seismic Wave Dynamics Basis[M]. Beijing:Petroleum Industry Press:1−2 (in Chinese).
童平. 2012. 地震层析成像方法及其应用研究[D]. 北京:清华大学:187−189. Tong P. 2012. Seismic Tomography Methods and Their Applications[D]. Beijing:Tsinghua University:187−189 (in Chinese).
王椿镛,段永红,吴庆举,王志铄. 2016. 华北强烈地震深部构造环境的探测与研究[J]. 地震学报,38(4):511–549. Wang C Y,Duan Y H,Wu Q J,Wang Z S. 2016. Exploration on the deep tectonic environment of strong earthquakes in North China and relevant research findings[J]. Acta Seismologica Sinica,38(4):511–549 (in Chinese).
殷文,印兴耀,吴国忱,梁锴. 2006. 高精度频率域弹性波方程有限差分方法及波场模拟[J]. 地球物理学报,49(2):561–568. Yin W,Yin X Y,Wu G C,Liang K. 2006. The method of finite difference of high precision elastic wave equations in the frequency domain and wave-field simulation[J]. Chinese Journal of Geophysics,49(2):561–568 (in Chinese).
张霖斌,姚振兴. 2000. 层状介质的声波波动方程反演[J]. 地球物理学进展,15(2):22–29. Zhang L B,Yao Z X. 2000. Wavform inversion of acoustic data in layered media[J]. Progress in Geophysics,15(2):22–29 (in Chinese).
张文生,罗嘉,滕吉文. 2015. 频率多尺度全波形速度反演[J]. 地球物理学报,58(1):216–228. doi: 10.6038/cjg20150119 Zhang W S,Luo J,Teng J W. 2015. Frequency multiscale full-waveform velocity inversion[J]. Chinese Journal of Geophysics,58(1):216–228 (in Chinese).
张文生,张丽娜. 2020. 基于有限元方法的频率域弹性波全波形反演[J]. 数值计算与计算机应用,41(4):315–336. Zhang W S,Zhang L N. 2020. Elastic wave full-waveform inversion based on the finite-element method in the frequency domain[J]. Journal on Numerical Methods and Computer Applications,41(4):315–336 (in Chinese).
祝贺君,刘沁雅,杨继东. 2023. 地震学全波形反演进展[J]. 地球与行星物理论评(中英文),54(3):287–317. Zhu H J,Liu Q Y,Yang J D. 2023. Recent progress on full waveform inversion[J]. Reviews of Geophysics and Planetary Phy sics,54(3):287–317 (in Chinese).
Ben-Hadj-Ali H,Operto S,Virieux J. 2011. An efficient frequency-domain full waveform inversion method using simultaneous encoded sources[J]. Geophysics,76(4):R109–R124. doi: 10.1190/1.3581357
Bornstein G,Biescas B,Sallarès V,Mojica J F. 2013. Direct temperature and salinity acoustic full waveform inversion[J]. Geophys Res Lett,40(16):4344–4348. doi: 10.1002/grl.50844
Brossier R,Operto S,Virieux J. 2009. Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[J]. Geophysics,74(6):WCC105–WCC118. doi: 10.1190/1.3215771
Bunks C,Saleck F M,Zaleski S,Chavent G. 1995. Multiscale seismic waveform inversion[J]. Geophysics,60(5):1457–1473. doi: 10.1190/1.1443880
Geng Y,Pan W Y,Innanen K A. 2018. Frequency-domain full-waveform inversion with non-linear descent directions[J]. Geophys J Int,213(2):739–756. doi: 10.1093/gji/ggy002
Hestholm S O,Ruud B O,Husebye E S. 1999. 3-D versus 2-D finite-difference seismic synthetics including real surface topography[J]. Phys Earth Planet Inter,113(1/2/3/4):339–354.
Jeong W,Lee H Y,Min D J. 2012. Full waveform inversion strategy for density in the frequency domain[J]. Geophys J Int,188(3):1221–1242. doi: 10.1111/j.1365-246X.2011.05314.x
Komatitsch D,Tromp J. 2003. A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation[J]. Geophys J Int,154(1):146–153. doi: 10.1046/j.1365-246X.2003.01950.x
Komatitsch D,Tsuboi S,Tromp J. 2005. The spectral-element method in seismology[G]//Seismic Earth:Array Analysis of Broadband Seismograms,Volume 157. Washington:American Geophysical Union:205.
Lang C,Yang D H. 2017. A nearly analytic discrete method for solving the acoustic-wave equations in the frequency domain[J]. Geophysics,82(1):T43–T57. doi: 10.1190/geo2016-0248.1
Lang C,Li Q S,Zhou Y J,He X J,Han R B. 2020. A high-precision low-dispersive nearly analytic difference method with its application in frequency-domain seismic waveform inversion[J]. Explor Geophys,51(3):355–377. doi: 10.1080/08123985.2019.1699786
Li A M,Liu H,Yuan Y X,Hu T,Guo X B. 2018. Modeling of frequency-domain elastic-wave equation with a general optimal scheme[J]. J Appl Geophys,159:1–15. doi: 10.1016/j.jappgeo.2018.07.014
Li J S,Yang D H,Liu F Q. 2013. An efficient reverse time migration method using local nearly analytic discrete operator[J]. Geophysics,78(1):S15–S23. doi: 10.1190/geo2012-0247.1
Liu S L,Yang D H,Ma J. 2017. A modified symplectic PRK scheme for seismic wave modeling[J]. Comput Geosci,99:28–36. doi: 10.1016/j.cageo.2016.11.001
Moczo P,Robertsson J O A,Eisner L. 2007. The finite-difference time-domain method for modeling of seismic wave propagation[J]. Adv Geophys,48:421–516.
Operto S,Virieux J,Amestoy P,L’Excellent J Y,Giraud L,Ali H B H. 2007. 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver:A feasibility study[J]. Geophysics,72(5):SM195–SM211. doi: 10.1190/1.2759835
Pratt R G,Worthington M H. 1988. The application of diffraction tomography to cross-hole seismic data[J]. Geophysics,53(10):1284–1294. doi: 10.1190/1.1442406
Pratt R G. 1990. Frequency-domain elastic wave modeling by finite differences:A tool for crosshole seismic imaging[J]. Geophysics,55(5):626–632. doi: 10.1190/1.1442874
Pratt R G,Shin C,Hick G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J]. Geophys J Int,133(2):341–362. doi: 10.1046/j.1365-246X.1998.00498.x
Pratt R G. 1999. Seismic waveform inversion in the frequency domain,part 1:Theory and verification in a physical scale model[J]. Geophysics,64(3):888–901. doi: 10.1190/1.1444597
Prieux V,Brossier R,Operto S,Virieux J. 2013. Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field,part 2:Imaging compressive-wave and shear-wave velocities[J]. Geophys J Int,194(3):1665–1681. doi: 10.1093/gji/ggt178
Pyun S,Son W,Shin C. 2009. Frequency-domain waveform inversion using an l1-norm objective function[J]. Explor Geophys,40(2):227–232. doi: 10.1071/EG08103
Qu Y M,Li J L,Li Z C,Huang J P. 2018. An elastic full-waveform inversion based on wave-mode separation[J]. Explor Geophys,49(4):530–552. doi: 10.1071/EG16158
Ren Z M,Liu Y. 2016. A hierarchical elastic full-waveform inversion scheme based on wavefield separation and the multistep-length approach[J]. Geophysics,81(3):R99–R123. doi: 10.1190/geo2015-0431.1
Saad Y. 2003. Iterative Methods for Sparse Linear Systems[M]. Philadelphia:Society for Industrial and Applied Mathematics:157−214.
Shin C,Yoon K,Marfurt K J,Park K,Yang D,Lim H Y,Chung S,Shin S. 2001. Efficient calculation of a partial-derivative wavefield using reciprocity for seismic imaging and inversion[J]. Geophysics,66(6):1856–1863. doi: 10.1190/1.1487129
Sirgue L,Pratt R G. 2004. Efficient waveform inversion and imaging:A strategy for selecting temporal frequencies[J]. Geophy sics,69(1):231–248. doi: 10.1190/1.1649391
Tape C,Liu Q Y,Tromp J. 2007. Finite‐frequency tomography using adjoint methods-methodology and examples using membrane surface waves[J]. Geophys J Int,168(3):1105–1129. doi: 10.1111/j.1365-246X.2006.03191.x
Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation[J]. Geophysics,49(8):1259–1266. doi: 10.1190/1.1441754
Tessmer E,Kosloff D,Behle A. 1992. Elastic wave propagation simulation in the presence of surface topography[J]. Geophys J Int,108(2):621–632. doi: 10.1111/j.1365-246X.1992.tb04641.x
Tong P,Yang D H,Hua B L. 2011. High accuracy wave simulation-revised derivation,numerical analysis and testing of a nearly analytic integration discrete method for solving acoustic wave equation[J]. Int J Solids Struct,48(1):56–70. doi: 10.1016/j.ijsolstr.2010.09.003
Tong P,Chen C W,Komatitsch D,Basini P,Liu Q Y. 2014. High-resolution seismic array imaging based on an SEM-FK hybrid method[J]. Geophys J Int,197(1):369–395. doi: 10.1093/gji/ggt508
Tromp J,Tape C,Liu Q Y. 2005. Seismic tomography,adjoint methods,time reversal and banana-doughnut kernels[J]. Geophys J Int,160(1):195–216.
Virieux J,Operto S. 2009. An overview of full-waveform inversion in exploration geophysics[J]. Geophysics,74(6):WCC1–WCC26. doi: 10.1190/1.3238367
Virieux J,Asnaashari A,Brossier R,Métivier L,Ribodetti A,Zhou W. 2017. An introduction to full waveform inversion[G]//Encyclopedia of Exploration Geophysics. Society of Exploration Geophysicists:R1-1−R1-40.
Wang Y H,Rao Y. 2006. Crosshole seismic waveform tomography,I:Strategy for real data application[J]. Geophys J Int,166(3):1224–1236. doi: 10.1111/j.1365-246X.2006.03030.x
Wang Y W,Dong L G,Liu Y Z,Yang J Z. 2016. 2D frequency-domain elastic full-waveform inversion using the block-diagonal pseudo-Hessian approximation[J]. Geophysics,81(5):R247–R259. doi: 10.1190/geo2015-0678.1
Yang D H,Teng J W,Zhang Z J,Liu E R. 2003. A nearly analytic discrete method for acoustic and elastic wave equations in anisotropic media[J]. Bull Seismol Soc Am,93(2):882–890. doi: 10.1785/0120020125
Yang D H,Lu M,Wu R S,Peng J M. 2004. An optimal nearly analytic discrete method for 2D acoustic and elastic wave equations[J]. Bull Seismol Soc Am,94(5):1982–1992. doi: 10.1785/012003155
Yang D H,Liu E R,Song G J,Wang N. 2009. Elastic wave modelling method based on the displacement-velocity fields:An improving nearly analytic discrete approximation[J]. J Seismol,13(2):209–217. doi: 10.1007/s10950-008-9122-2
Yang D H,Song G J,Zhang J H. 2010. A modified NAD algorithm with minimum numerical dispersion for simulation of anisotropic wave propagation[J]. J Seism Explor,19(1):21–42.
Yang J D,Zhu H J,Li X Y,Ren L,Zhang S. 2020. Estimating P wave velocity and attenuation structures using full waveform inversion based on a time domain complex‐valued viscoacoustic wave equation:The method[J]. J Geophys Res: Solid Earth,125(6):e2019JB019129. doi: 10.1029/2019JB019129
Youn O K,Zhou H W. 2001. Depth imaging with multiples[J]. Geophysics,66(1):246–255. doi: 10.1190/1.1444901
Zhao Z C,Chen J Y,Liu X B. 2020. Frequency-domain finite-difference elastic wave modeling in the presence of surface topography[J]. Pure Appl Geophys,177(6):2821–2839. doi: 10.1007/s00024-019-02402-1
-
本文视频
其他相关附件
-
PDF格式
2023-0038Appendix 点击下载(310KB)
-
计量
- 文章访问数: 228
- HTML全文浏览量: 60
- PDF下载量: 63
