2020, Vol. 63
2. 同济大学海洋与地球科学学院, 上海 200092;
3. 中交上海航道勘察设计研究院有限公司, 上海 200120
2. School of Ocean and Earth Science, Tongji University, Shanghai 200092, China;
3. Shanghai Waterway Engineering Co., Ltd., Shanghai 200120, China
全波形反演(Full Waveform Inversion, FWI)利用地震波波形信息揭示地下构造特征、获取储层物性，是研究地球内部或区域构造及油气资源勘探的常用手段.Tarantola(1984)最初提出了基于广义最小二乘理论的时域波形反演方法，此后全波形反演被推广至频域(Pratt, 1990, 1998; Pratt and Worthington 1990).目前，许多研究都集中于进一步开发、优化反演算法以提高效率和精度(Neumaier and Arnold, 1998; Pratt et al., 2002; Hu et al. 2011).
反问题和正问题密切联系.由于FWI依托于正演模拟获取合成地震记录，针对特定反演目标(如含起伏地表模型)，合适的正演方法能极大提高反演精度.有限差分法(Finite-Difference Method, FDM)以其易实现、计算高效等特点被广泛用于地震波模拟(Yerneni et al., 2002; Pei et al., 2012).但其采用规则网格剖分，会导致模型界面阶梯状离散而产生数值散射(Xue et al., 2014)，在界面崎岖不平或速度横向变化复杂时，影响尤为显著.针对此问题，交错网格(Madariaga, 1976; Levander, 1988)、可变网格(Moczo, 1989; Jastram and Tessmer, 1994; Huang and Dong, 2009a)和旋转交错网格(Saenger and Bohlen, 2004)等一系列算法被先后提出，但要在起伏地表下精确地实现自由边界条件仍是个复杂的问题.谱元法(Spectral Element Method, SEM)最初是在计算流体力学(Patera, 1984)中引入的.作为一种混合方法，其兼具有限元法可模拟任何复杂模型的韧性及伪谱法的高精度.Priolo和Seriani(1991)通过Chebyshev插值的1D声波谱元法实验证明，高阶插值函数可减少模拟所需的网格数量，提高了模拟效率.因此，采用非结构化网格以拟合起伏界面的谱元法，是模拟起伏自由地表波传播更好的选择.随后，勒让德多项式插值的谱元法也被成功应用于二维、三维地震波模拟及全球地震模拟当中(Komatitsch et al., 1999; Komatitsch and Tromp, 2002).
Tromp等(2005)在伴随状态法的框架下推导了目标函数关于密度、弹性参数梯度的表达式.在此基础上，谱元法被发展到2D/3D偏移成像、反演等领域(Tape et al., 2007, 2009; Fichtner et al., 2009; Luo et al., 2012, 2013; Borisov et al., 2017).而梯度的求取也可由散射积分法通过计算Fréchet导数直接得到，并能精确地实施预条件，但此时核函数的存储又会消耗大量内存，因此其很少在FWI中得到应用.Liu等(2015a)提出了改进的散射积分算法，实现在不存储核函数或直接使用Hessian矩阵下，计算一阶和二阶下降方向，目前已被成功地应用于全波形反演、逆时偏移和走时层析成像之中(Liu et al., 2015a; Yang et al. 2016a, b; Liu et al., 2018; Sun et al., 2017).
本文首先对比了谱元法和传统有限差分法针对起伏自由地表模型的正演效果，在此基础上，比较了基于两种正演引擎的FWI结果，研究了正演差异对反演效果的影响，并利用改进的散射积分算法优点，在频域实现了一种归一化的目标函数.最后，在理论分析稀疏观测系统对FWI的影响后，将基于谱元法的全波形反演方法应用于实际海底地震仪(Ocean Bottom Seismometer, OBS)数据的处理之中，获得了更多细节信息.
1 方法 1.1 正问题相比于传统有限差分法采用的规则网格，谱元法非结构化的网格剖分能更好地贴合起伏界面，减少数值散射(Fichtner, 2011).
以标量势ϕ表达的常密度声波方程(Peter et al., 2011)可表示为
|
(1) |
其中，x=(x, z)为空间位置，t表示时间，v表示速度，f表示震源，
|
(2) |
谱元法在方程(1)两边同时乘以任意一测试函数，并在整个模型空间Ω内积分得到波动方程的弱解形式.求解弱解形式可总结为三个步骤：(1)通过区域分解将模型空间分解为若干单元Ωe.二维情况下Ωe为不规则的四边形.再将各个单元映射到参考空间Λ=[-1, 1]×[-1, 1]上.参考空间Λ上的坐标系可表示为ξ=(ζ, η)；(2)通过n阶拉格朗日插值法，获得参考空间点上的波场值.波场φ(ξ, t)可由n+1个Gauss-Lobatto-Legende(GLL)采样点插值获得
|
(3) |
其中，lα和lβ为拉格朗日插值基函数，ξα和ηβ为GLL插值点，ϕ(ξα, ηβ, t)代表插值点(ξα, ηβ)上t时刻的波场值；(3)通过GLL求积法则将弱解方程中的积分转化为插值点上的累加求和.最终，弱解形式可写作关于时间的二阶微分方程：
|
(4) |
其中，M表示对角质量矩阵，K表示刚度矩阵，F为震源项.和最初的形式相比，方程(4)具有显式可解性.
为对比谱元法和传统有限差分法的模拟精度，本文设计一个具有起伏地表的两层模型(图 1).该模型以高斯函数刻画了一个高500 m的隆起和一个深300 m的凹陷，并包含一个具有300 m隆起的水平反射层.模型横宽7.5 km，纵深1.5 km.谱元法模拟中，非结构化网格水平方向82个元，垂直方向17个元，如图 2所示，元内采用4阶拉格朗日多项式插值.有限差分法以大小301×61，间隔25 m的规则网格离散模型，并使用10阶空间差分格式.炮点位于水平方向4.5 km，地表以下25 m处(由声波自由边界条件，地表处压力为0).
|
图 1 2D起伏地表模型 Fig. 1 The simple 2D model including topography |
图 3a和3b展示了相应的模拟结果.由于规则网格剖分导致起伏界面被阶梯状离散而形成虚假散射点，在有限差分法模拟结果中存在多余的散射波(红框).为减少数值散射，本文对空间间隔为5 m的细网格有限差分法进行测试.此时，模型被离散为1501×301个网格点，可以看到精细剖分差分法模拟记录与谱元法结果更为相似(图 4).图 5、6展示了从归一化的地震记录上抽取的水平1.25 km和5.375 km处的道记录，分别对应于地表和反射界面的隆起位置.结果表明，相比于谱元法，粗网格有限差分法(黑线)既存在数值散射(黄圈)又有相位上的差异(蓝色箭头).而细网格有限差分法(绿线)与谱元法结果吻合则较好.因此，对于传统有限差分法，若剖分过于稀疏，会导致模型界面离散不够准确，其形成的阶梯状结构将产生不必要的散射波和相位错误.
|
图 3 SEM(a)和粗网格FDM(b)模拟记录对比 Fig. 3 A comparison between simulation records for SEM (a) and FDM at 25 m grid spacing (b) |
|
图 4 SEM(a)和细网格FDM(b)模拟记录对比 Fig. 4 A comparison between simulation records for SEM (a) and FDM at 5 m grid spacing (b) |
|
图 5 x=1.25 km处道记录(左)和FDM与SEM记录之差(右) Fig. 5 Seismograms at x=1.25 km (left) and trace difference between FDM and SEM (right) |
|
图 6 x=5.375 km处道记录(左)和FDM与SEM记录之差(右) Fig. 6 Seismograms at x=5.375 km (left) and trace differences between FDM and SEM (right) |
表 1给出在不同模拟参数下，两种模拟方法的内存占用和计算耗时.虽然相比于粗网格差分法，谱元法计算时间是其4倍，内存占用几乎是其3倍，但是细网格差分法的内存占用几乎是谱元法的19倍，并且计算时间更长.有限元类型的方法通常被认为要求更多的内存和更高的计算量，该实验表明，在考虑起伏地表模拟精度的前提下，和有限差分法相比，谱元法实际上并非一种“昂贵”的方法.谱元法采用非结构化网格剖分，使得模拟更为精确，并避免产生人为散射波；这对严重依赖正演的FWI来说将起到非常重要的作用.
|
表 1 SEM与FDM内存消耗和计算耗时对比 Table 1 Comparisons of memory consumptions and CPU time between SEM and FDM |
由起伏界面不准确离散得到的阶梯状结构不仅会产生多余的散射波，同时也会引起相位移动.为分析正演效果对FWI的影响，本文对比了基于两种正演引擎的波形反演结果.
全波形反演旨在寻求合适的物理模型，使合成记录与观测记录达到最佳匹配.鉴于两种记录可能存在不同的能量等级，本文采用归一化目标函数(刘宝龙，2014)：
|
(5) |
其中，v表示模型参数，dc表示合成记录，dcref表示合成记录的参考道，d0表示观测记录，d0ref表示观测记录的参考道.若在子波非时变空变的假设下，该目标函数也可去除未知子波的影响.
经推导得到归一化目标函数的梯度：
|
(6) |
其中，
根据Woodward(1992)，声波方程的Fréchet偏导数可以表示为
|
(7) |
其中，g代表检波器的位置，s代表炮点位置，r代表地下任意空间点的位置，v0为背景速度，ω为角频率，G(g, r, ω)为检波点激发的格林函数，p(r, s, ω)为由炮点s激发的背景波场.若将g用gref代替之后，同理可得到参考道核函数.由此归一化数据的核函数K′可以表示为
|
(8) |
最终，目标函数(5)的梯度可以表示为
|
(9) |
本文FWI中采用改进的散射积分算法实现预条件最速下降方向的计算，即将以上矩阵-向量乘表达的梯度计算转化为具有明确物理含义的向量-标量乘的累加运算(式(10))
|
(10) |
避免了矩阵的显式生成与存储，既保证了Born正演算子和偏移算子的伴随性(互为共轭转置)，又大大降低了内存占用量，同时可以方便地实施预条件(式11).细节请参考Liu等(2015a).
|
(11) |
本文首先以图 1所示起伏地表模型进行基于两种正演引擎的波形反演对比实验.鉴于非结构化网格能更有效地贴合起伏界面，故以谱元法模拟记录作为波形反演的“观测数据”.两种方法的“合成记录”之间存在能量等级差异也印证使用归一化目标函数的必要性.反演的初始模型(图 7)为真实模型(图 1)高斯光滑结果.两种方法均在地表观测，在50~7450 m范围内均匀布设了149个炮点，间隔50 m，埋深25 m.每炮包含149个和炮点重合的检波点，并选择每炮零炮检距检波道作为当前炮道集的参考道.SEM正演的非结构化网格水平方向82个元，垂直方向17个元.FDM采用5 m细网格剖分.两种正演引擎的观测记录采样间隔均为0.4 ms.
|
图 7 简单模型实验的初始模型 Fig. 7 The starting model used in FWI test Ⅰ |
反演过程中将1~20 Hz的频率域数据不等间隔分为4个频率组.每个频率组最多迭代30次.值得一提的是，基于谱元法的全波形反演(SEM-FWI)是在地下反射界面位置未知的假设下开始的，故非结构化的网格仅贴合起伏地表.图 8为两种反演方法的第一代梯度.可见基于有限差分法的全波形反演(FDM-FWI)梯度在地表起伏处存在异常，投影到反演结果上(图 9)表现为界面起伏位置处速度不准确.水平位置1.25 km和2.25 km处的速度垂直抽线结果(图 10)也表明，FDM-FWI结果在起伏界面处存在抖动，而SEM-FWI则更稳定.两者在平地表处的结果基本相似.在效率方面，SEM-FWI和FDM-FWI所用时间比与其正演所花费的时间比一致，为1:19.6.
|
图 8 FDM-FWI(a)和SEM-FWI(b)的第一代梯度对比 Fig. 8 A comparison between first iteration gradients of FDM- (a) and SEM-based (b) FWIs |
|
图 9 FDM-FWI(a)和SEM-FWI(b)结果对比 Fig. 9 A comparison between the results of FDM- (a) and SEM-based (b) FWIs |
|
图 10 x=1.25 km(a)和x=2.25 km(b)处速度垂直抽线结果.黑线为真实模型速度，绿线为初始模型速度，红线表示SEM-FWI结果，蓝线表示FDM-FWI结果 Fig. 10 Vertical profiles at x=1.25 km (a) and at x=2.25 km (b). The black lines denote the velocity of the true model, the green lines represent the velocity of the starting model, the red lines indicate the velocity of the SEM-based FWI result, and the blue lines denote the velocity of the FDM-based FWI result. |
为验证本反演方法的可行性，进行含起伏地表的复杂模型(图 11)实验.该模型在Overthrust模型地表增加一个500 m的高斯隆起和300 m的凹陷.模型横宽20 km，纵深5 km.非结构化网格(图 12)仅贴合地表，水平方向223个元，垂直方向53个元，元内采用4阶多项式插值.在0~20 km之间均匀布设201个炮点，横向间隔100 m，埋深25 m.每炮包含201个与炮点重合的检波点.初始模型(图 13)为真实模型高斯光滑结果.反演过程中将2~20 Hz的频率域数据分为9个频率组.每组最大迭代次数30次.为检验模拟效果，同时进行基于有限差分法的全波形反演实验，共801×207个网格，间隔25 m，使用10阶空间差分格式.
|
图 11 Overthrust模型 Fig. 11 The true Overthrust model used in FWI test Ⅱ |
|
图 12 SEM非结构化网格 Fig. 12 The unstructured meshes used in SEM-based FWI |
|
图 13 Overthrust实验初始模型 Fig. 13 The starting model used in FWI test Ⅱ |
图 14对比了两种方法的反演结果，可以看出宏观的信噪比差异，图(b)比图(a)更加“干净”.局部上也可以发现，SEM-FWI在浅层更加稳定，FDM-FWI则存在错误的结构，如黑框所示；在中深层，前者空间结构也更为清晰，如蓝框所示.从水平抽线结果(图 15a和15b)也可见，FDM-FWI结果(蓝线)存在明显的抖动，而SEM-FWI结果(红线)更稳定；垂直抽线结果也表明，SEM-FWI结果更接近真实速度、具有更高分辨率，尤其是在3 km到4 km之间.
|
图 14 FDM-FWI(a)和SEM-FWI(b)结果对比 Fig. 14 A comparison between the results of FDM- (a) and SEM-based (b) FWIs |
|
图 15 z=200 m处(a)和z=2 km处(b)水平抽线；x=5 km处(c)和x=13.5 km处(d)垂直抽线.黑线为真实模型速度，绿线为初始模型速度，红线表示SEM-FWI结果，蓝线表示FDM-FWI结果 Fig. 15 Depth profiles 200 m below the surface (a), horizontal profiles at z=2 km (b), vertical profiles at x=5 km (c) and vertical profiles at x=13.5 km (d) |
OBS采集以其长偏移距、宽方位等优点被广泛用于地球深层构造研究和油气藏勘探.走时层析是处理OBS数据常用的方法(Katzman et al., 1994; Tinivella and Accaino, 2000; Gherasim et al., 2010)，但其忽略了丰富的波形信息.目前，利用OBS数据进行全波形反演仍处于理论研究阶段(Etienne et al., 2001; Sourbier et al., 2007; Kamei et al., 2013).本文将基于谱元法的全波形反演应用于西沙海槽OBS数据集以获取更多地学信息.
2.3.1 观测系统稀疏性影响理论实验OBS观测系统中检波点分布通常较为稀疏.为研究稀疏观测系统对反演的影响，在处理实际数据前应先进行理论测试.实验中采用了互易定理，将海底OBS点视为炮点，海面激发点视为检波点.图 16a为真实模型，横宽5.2 km，纵深1.5 km.谱元法非结构化网格水平剖分104个元，垂直剖分37个元.变形网格仅贴合海底，并在海面施加自由边界条件以模拟多次波.海水速度1.5 km·s-1.初始模型(图 16b)为真实模型光滑结果.
|
图 16 OBS理论实验采用的真实模型(a)和初始模型(b) Fig. 16 The true model (a) and the starting model (b) used in the OBS theoretical test |
本文中测试了两种观测系统.其一为稀疏观测系统.仅包含5个OBS，间隔1 km；其二为密集观测系统，共包含121个海底OBS，间隔40 m.图 17展示了反演结果.相比于密集观测系统，稀疏观测系统虽不能消除检波点间的采集脚印，并产生一系列规则抖动，但仍可以获得海底以下较准确的高波数结构.这说明稀疏观测系统下的实际资料全波形反演应用是可行的.
|
图 17 5个OBS(a)和121个OBS(b)反演结果 Fig. 17 FWI results with 5 OBSs (a) and 121 OBSs (b) respectively |
本文在南海西沙海槽的二维OBS数据上进行基于谱元法的全波形反演.实际观测系统如图 18所示.测线上包含14个间隔1 km的OBS和709个间隔50 m的震源.实验中仅利用了数据集中的声压分量.图 19a为OBSXS-1-2012经过炮检互易和基本处理后的共检波点道集.为提高信噪比，反演前须对数据进行一定的预处理：(1)球面扩散校正；(2)切除直达波前的噪声；(3)带通滤波去除40 Hz以上的能量；(4)FK滤波移除规则噪音.处理后的道集如图 19b所示.
|
图 18 西沙海槽观测系统.测线全长36 km，包含14个OBS(OBSXS-8和OBSXS-14丢失) Fig. 18 The geometry in the Xisha through. The line is 36 km in length with 14 OBSs (The OBSXS-8 and -14 were lost.) |
|
图 19 OBSXS-1-2012原始(a)和处理后(b)的共检波点道集 Fig. 19 Observed data (a) and processed data of OBS 1 (b) |
首先进行初至波走时层析，并将其结果作为全波形反演的初始模型(图 20)，其包含速度小于4.5 km·s-1的沉积层和速度小于8 km·s-1的上、下地壳(边界并不明显).谱元法非结构化网格如图 21所示，并且网格大小随着模型速度变化.实验中在海面施加自由边界条件，故不必从记录中去除多次波和鬼波.图 22展示了初始模型下于OBSXS-1-2012位置处模拟得到的共炮点道集.
|
图 20 由初至波走时层析获得的初始模型 Fig. 20 The starting model achieved from first-arrival traveltime tomography |
|
图 21 SEM采用的非结构化网格 Fig. 21 Meshes used in SEM modeling |
|
图 22 OBSXS-1-2012位置处基于初始模型的模拟记录 Fig. 22 Simulated data at the position of OBS 1 using SEM |
为了保证反演稳定性，在全波形反演中仅使用初至波波形信息.反演以谱元法为正演引擎合成模拟记录，并将观测数据和模拟数据分别进行带通滤波从而实施分频反演策略.数据被分为：4~8 Hz，8~12 Hz，12~15 Hz三个频带范围率组.每个频带数据的反演结果作为下一个频带数据的初始模型逐步反演.图 23展示了经加窗和滤波处理的观测数据及其相应的SEM-FWI目标函数下降曲线.
|
图 23 加窗、滤波观测数据(a)和各个频率组SEM-FWI目标函数下降曲线(b) Fig. 23 Windowed and filtered datasets (a) and misfits of each frequency group inversion (b) |
经过全波形反演之后，海底结构得到更新，沉积层和地壳中高波数异常变得明显(图 24).相比于初始模型，反演结果更加揭示了沉积层的横向非均质性.从模型更新量(图 25)中也可以看到这种效果.图 26展示了基于初始模型和反演结果的褶积数据残差，可以看到反演结果的数据残差有所减少(红圈)，尤其是折射波初至部分最为明显，证明了基于谱元法的全波形反演方法对OBS数据具有一定效果.值得一提的是，反演结果图 24中箭头所示的槽状构造与Shao等(2018)解释的渐新世早期形成的古河道(如图 27所示)位置上有很好的对应关系.
|
图 24 OBS数据反演结果 Fig. 24 The final FWI result for the field OBS dataset |
|
图 26 初始模型(a)和全波形反演结果(b)褶积残差(相比于观测数据的时间延迟由褶积造成) Fig. 26 Normalized data residuals based on the starting model (a) and the FWI result (b). The time shift compared with observed data is caused by convolution |
|
图 27 中南半岛至南海北部古河道(Shao et al., 2018).黄色短线表示OBS测线横穿该古河道 Fig. 27 The paleo channel from Indochina to North SCS (Shao et al., 2018) |
本文通过起伏地表模型正演实验说明，传统有限差分法对起伏界面规则剖分和阶梯状离散将导致数值散射和相移，而谱元法使用变形网格剖分能更有效地贴合起伏界面，因此具有较高的模拟精度.要获得与谱元法相似的模拟精度，有限差分法需要更精细的网格剖分，进而导致计算量上升和内存占用增大.故在对包含起伏界面模型进行模拟时，谱元法是一种更为合适的正演建模方法.
全波形反演的精度严重依赖于正演引擎.对含起伏地表模型的全波形反演实验表明，基于谱元法的波形反演能够重建更精细、稳定的地下速度结构.本文采用归一化的目标函数可以消除未知子波的影响，以及观测数据与模拟数据之间的能量差异.通过OBS观测系统的理论实验讨论了观测系统的稀疏性对反演的影响，虽然稀疏性会导致反演结果中存在一些规则抖动，但仍能获得海底以下的准确高波数结构.最后，本文将基于谱元法的全波形反演方法应用于实际OBS数据，并结合多尺度策略，揭示了西沙海槽地下深部更多的细节信息.
Borisov D, Modrak R, Gao F C, et al. 2017. 3D elastic full-waveform inversion of surface waves in the presence of irregular topography using an envelope-based misfit function. Geophysics, 83(1): R1-R11. |
Etienne G, Nicolétis L, Rakotoarisoa H. 2001. Comparison between different elastic imaging techniques on 2D OBS data.//71st Ann.Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1077-1080.
|
Fichtner A, Kennett B L N, Igel H, et al. 2009. Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods. Geophysical Journal International, 179(3): 1703-1725. DOI:10.1111/j.1365-246X.2009.04368.x |
Fichtner A. 2011. Full Seismic Waveform Modelling and Inversion. Berlin Heidelberg: Springer.
|
Gherasim M, Etgen J, Nolte B, et al. 2010. Anisotropic velocity model building using OBS node tomography at Atlantis field, Gulf of Mexico.//80th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 4359-4363.
|
Hu W, Abubakar A, Habashy T M, et al. 2011. Preconditioned non-linear conjugate gradient method for frequency domain full-waveform seismic inversion. Geophysical Prospecting, 59(3): 477-491. DOI:10.1111/j.1365-2478.2010.00938.x |
Huang C, Dong L G. 2009a. High-order finite-difference method in seismic wave simulation with variable grids and local time-steps. Chinese Journal of Geophysics (in Chinese), 52(1): 176-186. |
Huang C, Dong L G. 2009b. Staggered-grid high-order finite-difference method in elastic wave simulation with variable grids and local time-steps. Chinese Journal of Geophysics, 52(11): 2870-2878. |
Jastram C, Tessmer E. 1994. Elastic modelling on a grid with vertically varying spacing. Geophysical Prospecting, 42(4): 357-370. DOI:10.1111/j.1365-2478.1994.tb00215.x |
Kamei R, Pratt R G, Tsuji T. 2013. On acoustic waveform tomography of wide-angle OBS data-strategies for pre-conditioning and inversion. Geophysical Journal International, 194(2): 1250-1280. DOI:10.1093/gji/ggt165 |
Katzman R, Holbrook W S, Paull C K. 1994. Combined vertical-incidence and wide-angle seismic study of a gas hydrate zone. Blake Ridge. Journal of Geophysical Research:Solid Earth, 99(B9): 17975-17995. DOI:10.1029/94JB00662 |
Komatitsch D, Vilotte J P, Vai R, et al. 1999. The spectral element method for elastic wave equations-application to 2-D and 3-D seismic problems. International Journal for Numerical Methods in Engineering, 45(9): 1139-1164. DOI:10.1002/(SICI)1097-0207(19990730)45:9<1139::AID-NME617>3.0.CO;2-T |
Komatitsch D, Tromp J. 2002. Spectral-element simulations of global seismic wave propagation-Ⅰ. Validation. Geophysical Journal International, 149(2): 390-412. DOI:10.1046/j.1365-246X.2002.01653.x |
Levander A R. 1988. Fourth-order finite-difference P-SV seismograms. Geophysics, 53(11): 1425-1436. DOI:10.1190/1.1442422 |
Liu B L. 2014. Source independent full waveform inversion[Master's thesis] (in Chinese). Shanghai: Tongji University.
|
Liu Y Z, Yang J Z, Chi B X, et al. 2015a. An improved scattering-integral approach for frequency-domain full waveform inversion. Geophysical Journal International, 202(3): 1827-1842. DOI:10.1093/gji/ggv254 |
Liu Y Z, Wang G Y, Yang J Z, et al. 2015b. Multi-parameter full-waveform inversion for VTI media based on Born sensitivity kernels. Chinese Journal of Geophysics (in Chinese), 58(4): 1305-1316. |
Liu Y Z, Yang J Z, Zhou Q. 2018. An efficient nonlinear Fresnel volume tomography using an improved scattering-integral algorithm. Journal of Applied Geophysics, 159: 678-689. DOI:10.1016/j.jappgeo.2018.10.009 |
Luo Y, Tromp J, Denel B, et al. 2012. 3D elastic migration with topography based on spectral-element and adjoint methods.//82nd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1-5.
|
Luo Y, Tromp J, Denel B, et al. 2013. 3D coupled acoustic-elastic migration with topography and bathymetry based on spectral-element and adjoint methods. Geophysics, 78(4): S19-S202. |
Madariaga R. 1976. Dynamics of an expanding circular fault. Bulletin of the Seismological Society of America, 66(3): 639-666. |
Moczo P. 1989. Finite-difference technique for SH-waves in 2-D media using irregular grids-application to the seismic response problem. Geophysical Journal International, 99(2): 321-329. DOI:10.1111/j.1365-246X.1989.tb01691.x |
Neumaier A. 1998. Solving ill-conditioned and singular linear systems:a tutorial on regularization. SIAM Review, 40(3): 636-666. DOI:10.1137/S0036144597321909 |
Patera A T. 1984. A spectral element method for fluid dynamics:laminar flow in a channel expansion. Journal of Computational Physics, 54(3): 468-488. DOI:10.1016/0021-9991(84)90128-1 |
Pei Z L, Fu L Y, Sun W J, et al. 2012. Anisotropic finite-difference algorithm for modeling elastic wave propagation in fractured coalbeds. Geophysics, 77(1): C13-C26. |
Peter D, Komatitsch D, Luo Y, et al. 2011. Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes. Geophysical Journal International, 186(2): 721-739. DOI:10.1111/j.1365-246X.2011.05044.x |
Pratt R G, Worthington M H. 1990. Inverse theory applied to multi-source cross-hole tomography. Geophysical Prospecting, 38(3): 287-310. DOI:10.1111/j.1365-2478.1990.tb01846.x |
Pratt R G. 1990. Frequency-domain elastic wave modeling by finite differences:a tool for crosshole seismic imaging. Geophysics, 55(5): 626-632. DOI:10.1190/1.1442874 |
Pratt R G. 1998. Seismic waveform inversion in the frequency domain; Part 1, Theory and verification in a physical scale model. Geophysics, 64(3): 888-901. |
Pratt R G, Shin C, Hicks G J. 2002. Gauss-newton and full newton methods in frequency-space seismic waveform inversion. Geophysical Journal International, 133(2): 341-362. |
Priolo E, Seriani G. 1991. A numerical investigation of Chebyshev spectral element method for acoustic wave propagation.//Proceedings of the 13th IMACS World Congress on Computational and Applied Mathematics. Dublin: Criterion Press, 551-556.
|
Saenger E H, Bohlen T. 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid. Geophysics, 69(2): 583-591. |
Shao L, Cui Y C, Stattegger K, et al. 2018. Drainage control of Eocene to Miocene sedimentary records in the southeastern margin of Eurasian Plate. GSA Bulletin, 131(3-4): 461-478. |
Sourbier F, Operto S, Virieux J, et al. 2007. A massively parallel frequency-domain full-waveform inversion algorithm for imaging acoustic media: application to a dense OBS data set.//77th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1893-1897.
|
Sun M A, Yang J Z, Dong L G, et al. 2017. Density reconstruction in multiparameter elastic full-waveform inversion. Journal of Geophysics and Engineering, 14(6): 1445-1462. DOI:10.1088/1742-2140/aa93b0 |
Tape C, Liu Q Y, Tromp J. 2007. Finite-frequency tomography using adjoint methods-Methodology and examples using membrane surface waves. Geophysical Journal of the Royal Astronomical Society, 168(3): 1105-1129. DOI:10.1111/j.1365-246X.2006.03191.x |
Tape C, Liu Q, Maggi A, et al. 2009. Adjoint tomography of the southern California crust. Science, 325(5943): 988-992. DOI:10.1126/science.1175298 |
Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266. DOI:10.1190/1.1441754 |
Tinivella U, Accaino F. 2000. Compressional velocity structure and Poisson's ratio in marine sediments with gas hydrate and free gas by inversion of reflected and refracted seismic data (South Shetland islands, Antarctica). Marine Geology, 164(1-2): 13-27. DOI:10.1016/S0025-3227(99)00123-1 |
Tromp J, Tape C, Liu Q Y. 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels. Geophysical Journal International, 160(1): 195-216. |
Woodward M J. 1992. Wave-equation tomography. Geophysics, 57(1): 15-26. |
Xue Z, Dong L G, Li X B, et al. 2014. Discontinuous galerkin finite-element method for elastic wave modeling including surface topography. Chinese Journal of Geophysics, 57(4): 1209-1223. DOI:10.6038/cjg20140418 |
Yang J Z, Liu Y Z, Dong L G. 2016a. Least-squares reverse time migration in the presence of density variations. Geophysics, 81(6): S497-S509. DOI:10.1190/geo2016-0075.1 |
Yang J Z, Liu Y Z, Dong L G. 2016b. Simultaneous estimation of velocity and density in acoustic multiparameter full-waveform inversion using an improved scattering-integral approach. Geophysics, 81(6): R399-R415. DOI:10.1190/geo2015-0707.1 |
Yerneni S, Dheeraj D, Chakraborty S, et al. 2002. Finite difference forward modeling for complex geological models.//82nd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1987-1990.
|
黄超, 董良国. 2009a. 可变网格与局部时间步长的高阶差分地震波数值模拟. 地球物理学报, 52(1): 176-186. |
黄超, 董良国. 2009b. 可变网格与局部时间步长的交错网格高阶差分弹性波模拟. 地球物理学报, 52(11): 2870-2878. DOI:10.3969/j.issn.0001-5733.2009.11.022 |
刘玉柱, 王光银, 杨积忠, 等. 2015b. 基于Born敏感核函数的VTI介质多参数全波形反演. 地球物理学报, 58(4): 1205-1316. DOI:10.6038/cjg20150418 |
薛昭, 董良国, 李晓波, 等. 2014. 起伏地表弹性波传播的间断Galerking有限元数值模拟方法. 地球物理学报, 57(4): 1209-1223. DOI:10.6038/cjg20140418 |
刘宝龙. 2014.不依赖子波的全波形反演方法研究及应用[硕士论文].上海: 同济大学.
|