2019,
Vol. 43
2. 海洋国家实验室海洋矿产资源评价与探测技术功能实验室, 山东青岛 266071;
3. 中国石化中原油田分公司物探研究院, 河南濮阳 457001
2. Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China;
3. Geophysical Exploration Research Institute of Zhongyuan Oilfield Company, SINOPEC, Puyang 457001, China
全波形反演(FWI)作为一种较为重要的速度建模方法, 其能充分利用地震数据中所包含的振幅、相位、走时、波形等信息对地下模型参数进行精细刻画[1-6]。在全球和区域性构造的反演问题以及勘探地震中已取得一定的应用效果[7-8], 其中稳定性和计算成本是全波形反演(FWI)的两个重要的研究方面[9-12]。在反演过程中, 给定的是精确或者近似精确的Hessian矩阵, 则不需要进行步长的计算[13], 但是很难得到精确的Hessian矩阵, 一般应用对角阵或者伪Hessian矩阵进行近似, 因此需要计算附加步长[14-16], 其中抛物插值思想求取步长是一种广泛使用且有效的方法[17], 其需要两步或者多步向前正演模拟, 计算量巨大。在保证地震波场正演模拟精度的同时, 如何提高正演计算效率是减少波形反演计算成本的关键所在。刘玉柱等[18]利用高斯束计算格林函数和正演波场, 减少了正演计算量, Berkhout[19]提出的同时激发震源技术, 将不同空间位置的多个震源按照随机线性编码方式激发组成超级炮。Krebs、张生强等[20-21]在全波形反演过程中使用了同时激发震源技术, 结合编码技术减少串扰噪音, 但是其本质上不能够消除串扰噪。Wang等[22]提出了利用编码的多震源方法直接求取步长, 其在应用中需要给定初始系数, 并且系数不同计算的步长也不一样。笔者在抛物插值的基础上, 将多震源混叠数据反演的思想引入步长计算中, 实现基于混叠数据的优化步长声波全波形反演方法。
1 原理 1.1 时间域全波形反演基本原理传统基于L2范数的全波形反演方法的目标泛函为
| $ E(m) = \sum\limits_{n = 1}^{{N_{\rm{s}}}} {{{\left\| {{u_{{\rm{cal}}}}\left( {m,{s_n}} \right) - {u_{{\rm{obs}}n}}} \right\|}^2}} . $ | (1) |
式中, uobs、ucal分别为由震源sn激发的观测数据和正演数据; m为模型参数; Ns为总炮数。
目标泛函关于模型的梯度项可以通过伴随状态法求得[23]:
| $ {\nabla _m}E\delta m = \left\langle {\hat u{\nabla _m}B\left( u \right)\delta m} \right\rangle . $ | (2) |
式中,
| $ B(m) u=S. $ | (3) |
以及
| $ \hat{B}(m) \hat{u}=S_{\mathrm{adj}}. $ | (4) |
式中, S和Sadj分别为震源和伴随震源; 为伴随正演算子。伴随震源Sadj可以表示为
| $ S_{\mathrm{adj}}=u_{\mathrm{cal}}-u_{\mathrm{obs}}. $ | (5) |
全波形反演通过对模型参数的迭代更新求得最优解, 其迭代公式可以表示为
| $ m_{k+1}=m_{k}-\alpha_{k} H_{k} \nabla_{m} E. $ | (6) |
式中, k为迭代次数; Hk为对梯度的修正项, 可采用共轭梯度法、L-BFGS方法等[21, 24-25], 本文中采用共轭梯度法; αk为步长。
1.2 抛物步长求取算法步长计算是保证反演效率和精度的关键因素。传统反演中多采用试探步长法, 本文中选取抛物插值法作为标准的方法进行测试(traditional parabolic interpolation method, TPI)分别选取两个试算步长α1和α2, 计算其对应的目标泛函值, 如图 1所示。为了更好地比较计算量和反演误差, 本文中采用两个试探步长, 则最优步长可以利用下式获得:
|
图 1 抛物拟合求取步长示意图[17] Fig.1 Diagram of parabolic fitting calculating step-length |
| $ {\alpha _{{\rm{best}}}} = \frac{1}{2}\frac{{\left( {{E_1} - {E_0}} \right)\alpha _2^2 - \left( {{E_2} - {E_0}} \right)\alpha _1^2}}{{\left( {{E_1} - {E_0}} \right){\alpha _2} - \left( {{E_2} - {E_0}} \right){\alpha _1}}}. $ | (7) |
式中, E0为第k次迭代的误差泛函; E1、E2为试探步长α1、α2的误差泛函, 将式(1)代入, 可以得到传统的步长计算公式:
| $ {a_{{\rm{opt}} - {\rm{trd}}}} = \frac{1}{2}\frac{{\left( {\sum\limits_{n = 1}^{{N_{\rm{s}}}} {{{\left\| {{u_{{\rm{cal}},n}}\left( {{m_1},{s_n}} \right) - {u_{{\rm{obs}},n}}} \right\|}^2}} - {E_k}} \right)\alpha _2^2 - \left( {\sum\limits_{n = 1}^{{N_{\rm{s}}}} {{{\left\| {{u_{{\rm{cal}},n}}\left( {{m_2},{s_n}} \right) - {u_{{\rm{obs}},n}}} \right\|}^2}} - {E_k}} \right)\alpha _1^2}}{{\left( {\sum\limits_{n = 1}^{{N_{\rm{s}}}} {{{\left\| {{u_{{\rm{cal}},n}}\left( {{m_1},{s_n}} \right) - {u_{{\rm{obs}},n}}} \right\|}^2}} - {E_0}} \right){\alpha _2} - \left( {\sum\limits_{n = 1}^{{N_{\rm{s}}}} {{{\left\| {{u_{{\rm{cal}},n}}\left( {{m_2},{s_n}} \right) - {u_{{\rm{obs}},n}}} \right\|}^2}} - {E_0}} \right){\alpha _1}}}. $ | (8) |
其中
| $ m_{1}=m_{k}-\alpha_{1} H_{k} \nabla_{m} E, $ |
| $ m_{2}=m_{k}-\alpha_{2} H_{k} \nabla_{m} E, $ |
| $ E_{k}=\sum\limits_{n=1}^{N_{\mathrm{s}}}\left\|u_{\mathrm{cal}, n}\left(m_{k}, s_{n}\right)-u_{\mathrm{obs}, n}\right\|^{2} $ |
如果严格按照抛物插值的条件则需要进行另外的试探步长计算, 大大占用了计算资源。当出现计算步长为负值或过大时, 本文中采取以下的计算策略:
| $ {a_{{\rm{opt - trd}}}} = \left\{ {\begin{array}{*{20}{l}} {{\alpha _{l\max }},{a_{{\rm{cal}},{\rm{opt}}}} > {\alpha _{l\max }};}\\ {{\alpha _1},{a_{{\rm{cal}},{\rm{opt}}}} < 0.} \end{array}} \right. $ | (9) |
式中, acal, opt为计算的步长; aopt-trd为采用的步长; almax为根据实际情况限定的最大步长。
1.3 多震源并行步长估计由于采用两个试探步长, 其正演需要巨大的计算量, 为了提高计算效率, 本文中提出多震源抛物插值算法(multi-sourceparabolic interpolation method, MSPI)求取参数修正步长。其选取3个试算步长α1、α2和α3, 采用多震源正演估算步长的误差泛函值, 可以大幅减少试探步长的正演炮数, 由于采用多震源激发, 一次步长计算时间相当于一炮正演的时间, 可以减少计算时间。基于多震源估计的抛物步长公式为
| $ \begin{array}{l}{a_{\mathrm{opt}-\mathrm{mult}}=} \\ {-\frac{1}{2} \frac{\left(\alpha_{2}^{2}-\alpha_{1}^{2}\right)\left(E_{m 3}-E_{m 1}\right)-\left(\alpha_{3}^{2}-\alpha_{1}^{2}\right)\left(E_{m 2}-E_{m 1}\right)}{\left(\alpha_{3}-\alpha_{1}\right)\left(E_{m 2}-E_{m 1}\right)-\left(\alpha_{2}-\alpha_{1}\right)\left(E_{m 3}-E_{m 1}\right)}}\end{array}. $ | (10) |
其中
| $ E_{m 1}=\left\|u_{\mathrm{cal}, n}\left(m_{1}, \sum\limits_{n=1}^{N_{\mathrm{s}}} s_{n}\right)-\sum\limits_{n=1}^{N_{\mathrm{s}}} u_{\mathrm{obs}, n}\right\|^{2}, $ |
| $ E_{m 2}=\left\|u_{\mathrm{cal}, n}\left(m_{2}, \sum\limits_{n=1}^{N_{\mathrm{s}}} s_{n}\right)-\sum\limits_{n=1}^{N_{\mathrm{s}}} u_{\mathrm{obs}, n}\right\|^{2}, $ |
| $ E_{m 3}=\left\|u_{\mathrm{cal}, n}\left(m_{3}, \sum\limits_{n=1}^{N_{\mathrm{s}}} s_{n}\right)-\sum\limits_{n=1}^{N_{\mathrm{s}}} u_{\mathrm{obs}, n}\right\|^{2}, $ |
| $ m_{1}=m_{k}-\alpha_{1} H_{k} \nabla_{m} E, $ |
| $ m_{2}=m_{k}-\alpha_{2} H_{k} \nabla_{m} E, $ |
| $ m_{3}=m_{k}-\alpha_{3} H_{k} \nabla_{m} E. $ |
由于MSPI采用3个试探步长(传统方法中, 步长为零算作一次试探步长), 由于实际误差以及目标泛函不匹配的影响, 选取的试探步长不准确, 尤其在初始模型准确时产生的误差较大, 导致收敛速度较慢, 对收敛结果产生影响。为了更加准确有效地估算步长, 本文中提出多震源多个步长优选方法(multi-source multi-step optimization, MSMO), 对MSPI方法进行改进。
由于采用多震源正演的算法, 其计算过程中采用并行计算, 这样只进行3个试探步长估计的MSPI方法, 在实际的计算中对计算资源有巨大的闲置, 所以本文中采用MSMO方法, 即选取Nu个试探步长进行计算(本文中在计算时选取Nu=10), 并从中选取最优的步长为迭代步长。
| $ \alpha_{\mathrm{opt-xu}}=\min \left\{E_{m 1}, E_{m 2}, \cdots, E_{m N_{\mathrm{u}}}\right\}. $ | (11) |
其中
| $ {E_{m1}} = {\left\| {{u_{{\rm{cal}},n}}\left( {{m_1},\sum\limits_{n = 1}^{{N_{\rm{s}}}} {{s_n}} } \right) - \sum\limits_{n = 1}^{{N_{\rm{s}}}} {{u_{{\rm{obs}},n}}} } \right\|^2}, $ |
| $ {E_{m2}} = {\left\| {{u_{{\rm{cal}},n}}\left( {{m_2},\sum\limits_{n = 1}^{{N_{\rm{s}}}} {{s_n}} } \right) - \sum\limits_{n = 1}^{{N_{\rm{s}}}} {{u_{{\rm{obs}},n}}} } \right\|^2}, $ |
| $ {E_{m{N_{\rm{u}}}}} = {\left\| {{u_{{\rm{cal}},n}}\left( {{m_{{N_{\rm{u}}}}},\sum\limits_{n = 1}^{{N_{\rm{s}}}} {{s_n}} } \right) - \sum\limits_{n = 1}^{{N_{\rm{s}}}} {{u_{{\rm{obs}},n}}} } \right\|^2}, $ |
| $ m_{1}=m_{k}-\alpha_{1} H_{k} \nabla_{m} E, $ |
| $ {m_2} = {m_k} - {\alpha _2}{H_k}{\nabla _m}E, $ |
| $ {m_{{N_{\rm{u}}}}} = {m_k} - {\alpha _{{N_{\rm{u}}}}}{H_k}{\nabla _m}E. $ |
由于计算了多个步长, 可以更加有效地估计步长, 特别是在前期步长波动较大的时候, MSPI插值方法存在较大误差, 采用MSMO方法可以更好地更新速度场, 从而在减少计算量同时, 保持了TPI较快的收敛速度和较高的反演精度, 从而提高反演计算效率。
在实际的步长估计中可以用以下的公式确定EmNu对应的步长:
| $ {\alpha _{m{N_{\rm{u}}},k}} = \left\{ {\begin{array}{*{20}{l}} {{\alpha _0},k = 1;}\\ {2{a_{{\rm{opt}} - {\rm{xu}},k - 1}} + \delta (k),k > 1.} \end{array}} \right. $ | (12) |
式中, αmNu, k为第k次迭代选取步长的范围的最大值; α0为反演中初始步长根据实际情况可选取的最大步长; αopt-xu, k-1为第k-1次迭代中选取的最优步长; δ(k)为0到α0/10的一个随机数, 防止后期步长更新陷入一个常数。选取的Nu个试探步长可以在允许的最小步长αmin到αmNu, k之间均匀选取, 也可以随机选取, 本文中在试算时采用均匀选取策略。
1.4 计算量及计算时间估算本文中方法可以大大提高计算效率, 为了估计3种方法的计算量, 设传统的全波形反演单炮一次正演所需时间记为T0, 总炮数为Ns, 则3种方法一次迭代的近似计算量分别为:传统方法为4NsT0; 多震源抛物插值为(2Ns+3)T0; 多震源优选为(2Ns+Nu)T0。
1.5 基于多震源估计的优化步长声波FWI流程从基于多震源估计的优化步长声波全波形反演方法流程图 2可知:基于地质与地球物理观测数据建立合适的背景速度模型, 应用有限差分进行正演得到炮记录。构建一个合适的目标泛函, 对每一炮计算每检波点位置的波场残差, 将波场残差反传播到模型空间, 得到反传播数据, 求取反传波场与正传波场在时间上的二阶导数关于时间变量的内积, 得到单炮的梯度。叠加所有炮求取的梯度值, 得到模型空间的全局梯度, 给定多个全局的试探步长(如将速度更新量控制在背景速度的1%), 以此更新速度得到一个伪更新速度模型, 进行多震源同步正演得到混叠记录, 对观测数据进行叠加得到超级炮数据, 求取多震源误差泛函, 然后利用误差泛函计算合适的步长, 进行速度更新, 直到迭代得到合适的速度场。
|
|
图 2 基于多震源多步长优选的全波形反演流程 Fig.2 Flow chart of FWI based on multi-source multi-step optimization |
为了验证方法的有效性, 本文中先应用简单模型(图 3(a))进行测试, 将上述3种步长估算方法进行模型测试, 纵横向尺寸为2.5 km×4.26 km, 每200 m做一次正演模拟, 震源子波主频为15 Hz, 采样时长为2.4 s。本文中使用时间2阶, 空间8阶有限差分正演模拟[26]及PML边界条件[27-29]。放炮方式采用地表全接收, 共40炮, 炮间距为100 m。采用三角平滑方式(0.5 km×1.5 km)得到反演的初始模型(图 3(b))所示, 本文中先采用时间域全波形算法进行测试, 反演过程中给定震源子波。
|
图 3 反演采用的速度模型 Fig.3 Velocity model for inversion |
3种步长估算方法的反演结果如图 4所示。3种方法反演结果与真实模型相比, 速度数值较为接近, 主要层位置及断层形态都得到了较好恢复。多震源估计步长的两种方法(MSPI、MSMO)相比于传统的方法(TPI), 反演结果相似, 说明其具有较好的反演精度, 估算的步长较为准确。为了更好地展示3种方法的差异性, 本文中给出了3种方法的收敛曲线, 如图 5所示。通过收敛曲线, 可以看出3种步长的估算方法都能较好收敛, 不同的估算方法在收敛速度和程度上只在开始阶段存在明显差异。在初始速度场比较准确的时候, 步长对速度场的更新影响较小, 基于多震源的步长估算方法(MSPI、MSMO)都能较好收敛, 与传统的方法(TPI)收敛曲线近似一致。
|
图 4 不同方法的反演结果 Fig.4 Inversion results by different methods |
|
图 5 不同全波形反演方法的数据误差 Fig.5 Data residuals of different FWIs |
本文中给出的两种多震源估计步长的方法在计算效率上有显著的提高, 其计算量如图 6所示(图中计算量表示不同方法需要计算量与一次单炮正演计算量的比值), 本文中采用10个节点并行计算。计算结果表明:传统全波形反演方法(TPI)求取的步长计算量是多震源估计步长方法(MSPI、MSMO)的近10倍, 总计算是多震源估计步长方法的近两倍。采用多震源估计的方法能大大减少步长计算量, 而MSMO的方法计算量稍大于MSPI的方法。两种多震源估计步长的方法, 在保持反演精度的同时, 可以减少用于步长的计算量, 显著提高反演的计算效率, 这样对于全波形的实用化具有一定意义。
|
图 6 不同反演方法的计算量比较 Fig.6 Computational costs for different FWIs |
在全波形反演中如果初始速度模型与真实速度模型相差较大, 而反演的速度模型又比较复杂情况下, 单尺度的全波形方法存在周波跳跃现象, 采用多尺度的全波形反演方法是解决这一问题比较有效地方法[30]。由于其采用分频的计算方法, 所以多尺度全波形反演方法的计算量相比于单尺度方法显著增加, 而采用多震源估算步长可以有效地减少求取步长计量, 其具有重要的意义。为验证本文中两种多震源估算步长方法的适用性, 将其应用到复杂的Marmousi模型(图 7(a)), 初始速度场为线性增加的梯度场(图 7(b)), 纵横向尺寸为1.5 km×4.58 km, 震源子波主频为12 Hz, 采样时长为4.0 s。采用全接收炮方式, 共44炮, 炮间距为100 m。为了能够保证照明均匀, 对模型右侧进行扩边处理(扩边部分未显示), 模型存在50 m的水层(图 7), 可见初始模型与真实模型相差较大。反演过程中假定震源子波已知, 水层速度已知, 采用基于维纳滤波的多尺度声波时间域全波形反演方法[30]。
|
图 7 多尺度反演采用的速度模型 Fig.7 Vvelocity model for multi-scale inversion |
本文中采用从2 Hz主频开始反演, 每迭代10次增加1 Hz的主频反演策略。第50次迭代(主频为2~6 Hz)3种步长估算方法的反演结果如图 8所示。传统的多尺度全波形方法(TPI)可以较好地对背景场进行反演, 复杂区域的大构造相对准确, 与真实值更加接近; 而多震源抛物拟合方法(MSPI)与TPI相比, 可以看出由于初始模型与真实模型相差太远, 反演陷入了局部极小值, 不能收敛到正确结果, 求取的步长不准确, 导致收敛速度慢, 只恢复出了浅部构造的大致轮廓, 而中深部及复杂区域与真实模型差别较大; 而多震源优选步长反演方法(MSMO)与TPI的反演结果相近, 说明在初始反演阶段MSMO的方法相比于MSPI求取的步长比较合理,也说明步长对于全波形反演具有重要意义。3种步长估算方法的多尺度全波形反演方法经过130次迭代的反演结果如图 9所示, 可以看出MSPI反演方法(图 9(b))在浅部有些区域陷入了局部极小值, 对于复杂区域的层位不能准确刻画, 这是由于前期收敛速度较慢导致反演背景场不准确, 使得后期反演容易陷入局部极值, 另外对深部的反演也比较弱。TPI(图 9(a))与MSMO(图 9(c))两种方法对于复杂模型具有更强的重建能力, 浅中部界面清晰, 深部构造也比较明显, 反演得到的速度更加接近真实值。MSMO方法在反演的精度上近似接近于TPI方法, 但是其用于步长的计算量大大减少, 显著提高了反演的计算效率。
|
图 8 Marmousi模型的50次迭代反演结果 Fig.8 Inversion profile of Marmousi velocity model after 50 iterations |
|
图 9 Marmousi模型的130次迭代反演结果 Fig.9 Inversion profile of Marmousi velocity model after 130 iterations |
为了更好地比较3种方法的差异性, 本文中给出了3种方法的收敛曲线, 如图 10所示。通过收敛曲线, 可以看出TPI传统的抛物拟合方法能较好地收敛到一个比较小的数值。MSPI方法的误差曲线也可以收敛但是收敛速度较慢, 并且最后收敛的数值相对于TPI较大, 这是由于步长估算不准确, 导致前期更新较慢, 背景场反演不准确, 导致后期陷入局部极值。MSMO方法也可以较好地收敛, 与TPI的误差收敛曲线相接近, 其略差于TPI, 这是由采用多震源估算步长时存在一些误差造成的。但是相比于传统的TPI方法, MSMO方法可以节省大量的步长计算量。本文中采用44炮计算, 则在计算步长方面, MSMO方法的步长计算量约为TPI的1/9, 如果炮数增加, 则计算效率的提高会更明显。
|
图 10 不同尺度反演方法的数据误差 Fig.10 Data residuals of different multi-scale inversion |
(1) 采用传统的抛物拟合的方法(TPI)可以有效地估算步长, 误差曲线可以很好地收敛, 并且反演精度较高, 但是要对试探步长更新的速度场进行两次正演, 计算量巨大。
(2) 采用多震源抛物插值方法(MSPI)在初始速度相对准确时可以较好地估算步长, 反演精度和收敛曲线与传统方法基本相同, 其在步长估计过程中采用多震源正演, 大大减少正演的计算量, 极大地提高了计算效率, 但是MSPI方法对于初始模型与真实模型相差较大的情况下适用性不好。
(3) 采用多震源多步长优选方法(MSMO)估算步长, 其对于简单和复杂模型都具有较好的适应性, 反演结果与TPI方法相接近, 具有较好的收敛性和稳定性, 但是相比于TPI方法, 其在步长估计过程中采用多震源正演, 并且采用并行技术只需要一次正演时间, 极大地提高了计算效率, 并且其对多尺度反演方法具有较好的适用性。
| [1] |
TARANTOLA A. Inversion of seismic reflection data in the acoustic approximation[J]. Geophysics, 1984, 49(8): 1259-1266. DOI:10.1190/1.1441754 |
| [2] |
VIRIEUX J, OPERTO S. An overview of full-waveform inversion in exploration geophysics[J]. Geophysics, 2009, 74(6): WCC1-WCC26. |
| [3] |
SHIN C, CHA Y H. Waveform inversion in the Laplace domain[J]. Geophysical Journal of the Royal Astronomical Society, 2008, 173(3): 922-931. DOI:10.1111/gji.2008.173.issue-3 |
| [4] |
黄建平, 崔超. 一种基于动态目标泛函的全波形反演多解性评估方法[J]. 中国石油大学学报(自然科学版), 2017, 41(6): 64-70. HUANG Jianping, CUI Chao. Multi-solution analysis of full waveform inversion based on dynamic misfit function[J]. Journal of China University of Petroleum (Edition of Natural Science), 2017, 41(6): 64-70. DOI:10.3969/j.issn.1673-5005.2017.06.007 |
| [5] |
杨积忠, 刘玉柱, 董良国. 变密度声波方程多参数全波形反演策略[J]. 地球物理学报, 2014, 57(2): 628-643. YANG Jizhong, LIU Yuzhu, DONG Liangguo. A multi-parameter full waveform inversion strategy for acoustic media with variable density[J]. Chinese Journal of Geophysics, 2014, 57(2): 628-643. |
| [6] |
MORA P. Nonlinear two-dimensional elastic inversion of multi-offset seismic data[J]. Geophysics, 1987, 52(9): 1211-1228. DOI:10.1190/1.1442384 |
| [7] |
FICHTNER A, BRIAN L N, HEINER I, et al. Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods[J]. Geophysical Journal International, 2009, 179(3): 1703-1725. DOI:10.1111/gji.2009.179.issue-3 |
| [8] |
TAPE C, LIU Q, MAGGI A, et al. Adjoint tomography of the Southern California crust[J]. Science, 2009, 325(5943): 988-992. DOI:10.1126/science.1175298 |
| [9] |
WANG H, SINGH S C, AUDEBERT F, et al. Inversion of seismic refraction and reflection data for building long-wavelength velocity models[J]. Geophysics, 2015, 80(2): R81-R93. |
| [10] |
BROSSIER R, OPERTO S, VIRIEUX J. Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[J]. Geophysics, 2009, 74(6): WCC105-WCC118. |
| [11] |
TARANTOLAA. A strategy for nonlinear elastic inversion of seismic reflection data[J]. Geophysics, 1986, 51(10): 1893-1903. DOI:10.1190/1.1442046 |
| [12] |
BUNKS C, SALECK F M, ZALESKI S, et al. Multiscale seismic waveform inversion[J]. Geophysics, 1995, 60(5): 1457-1473. |
| [13] |
FICHTNER A, TRAMPERT J. Hessian kernels of seismic data functionals based upon adjointtechniques[J]. Geophysical Journal International, 2011, 185(2): 775-798. DOI:10.1111/gji.2011.185.issue-2 |
| [14] |
CHOI Y, MIN D J, SHIN C. Frequency-domain elastic full waveform inversion using the new pseudo-Hessian matrix: experience of elastic Marmousi-2 synthetic data[J]. Bulletin of the Seismological Society of America, 2008, 98(5): 2402-2415. DOI:10.1785/0120070179 |
| [15] |
GAUTHIER O, VIRIEUX J, TARANTOLA A. Two-dimensional nonlinear inversion of seismic waveforms: numerical results[J]. Geophysics, 1986, 51(7): 1387-1403. DOI:10.1190/1.1442188 |
| [16] |
WRIGHT S, NOCEDAL J. Numerical optimization[J]. Springer Science, 1999, 35(67/68): 7. |
| [17] |
VIGH D, STARR E W, KAPOOR J. Developing earth models with full waveform inversion[J]. The Leading Edge, 2009, 28(4): 432-435. |
| [18] |
LIU Y, XIE C, YANG J. Gaussian beam first-arrival waveform inversion based on Born wavepath[J]. Chinese Journal of Geophysics, 2014, 57(5): 619-630. DOI:10.1002/cjg2.2014.57.issue-5 |
| [19] |
BERKHOUT A J G. Changing the mindset in seismic data acquisition[J]. The Leading Edge, 2008, 27(7): 924-938. DOI:10.1190/1.2954035 |
| [20] |
KREBS J R, ANDERSON J E, HINKLEY D, et al. Fast full-wavefield seismic inversion using encoded sources[J]. Geophysics, 2009, 74(6): WCC177-WCC188. |
| [21] |
张生强, 刘春成, 韩立国, 等. 基于L-BFGS算法和同时激发震源的频率多尺度全波形反演[J]. 吉林大学学报(地球科学版), 2013, 43(3): 1004-1012. ZHANG Shengqiang, LIU Chuncheng, HAN Liguo, et al. Frequency multi-scale full waveform inversion based on L-BFGS algorithm and simultaneous sources approach[J]. Journal of Jilin University(Earth Science Edition), 2013, 43(3): 1004-1012. |
| [22] |
WANG S, CHEN J. Step length of full-waveform inversion based on linear search method in multisource domain[C/OL]//2016 SEG Annual Meeting, Dallas, Texas, October 16-21, 2016: SEG Technical Program Expanded Abstracts, c2016: 1492-1495.[2017-03-10]. https: //library.seg.org/doi/10.1190/segam2016-13771939.1.
|
| [23] |
PLESSIX R E. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[J]. Geophysical Journal International, 2010, 167(2): 495-503. |
| [24] |
BECK A, TEBOULLE M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J]. Siam Journal on Imaging Sciences, 2009, 2(1): 183-202. DOI:10.1137/080716542 |
| [25] |
王义, 董良国. L-BFGS法时间域全波形反演中初始矩阵的选择方法[J]. 石油物探, 2014, 53(5): 545-555. WANG Yi, DONG Liangguo. Selection strategy of the initial matrix for L-BFGS method in time domain full waveform inversion[J]. Geophysical Prospecting for Petroleum, 2014, 53(5): 545-555. DOI:10.3969/j.issn.1000-1441.2014.05.007 |
| [26] |
ALAN R L. Fourth-order finite-difference P-SV seismograms[J]. Geophysics, 1988, 53(11): 425-1436. |
| [27] |
田坤, 黄建平, 李振春, 等. 多轴卷积完全匹配层吸收边界条件[J]. 石油地球物理勘探, 2014, 49(1): 143-152. TIAN Kun, HUANG Jianping, LI Zhenchun, et al. Multi-axial convolution perfectly matched layer absorption boundary condition[J]. Oil Geophysical Prospecting, 2014, 49(1): 143-152. |
| [28] |
黄建平, 杨宇, 李振春, 等. 基于M-PML边界的Lebedev网格起伏地表正演模拟方法及稳定性分析[J]. 中国石油大学学报(自然科学版), 2016, 40(4): 47-56. HUANG Jianping, YANG Yu, LI Zhenchun, et al. Lebedev grid finite difference modeling for irregular free surface and stability analysis based on M-PML boundary condition[J]. Journal of China University of Petroleum (Edition of Natural Science), 2016, 40(4): 47-56. DOI:10.3969/j.issn.1673-5005.2016.04.006 |
| [29] |
雍鹏, 黄建平, 李振春, 等. 优化的时空域等效交错网格有限差分正演模拟[J]. 中国石油大学学报(自然科学版), 2017, 41(6): 71-79. YONG Peng, HUANG Jianping, LI Zhenchun, et al. Forward modeling by optimized equivalent staggered-grid finite-difference method for time-space domain[J]. Journal of China University of Petroleum (Edition of Natural Science), 2017, 41(6): 71-79. DOI:10.3969/j.issn.1673-5005.2017.06.008 |
| [30] |
BOONYASIRIWAT C, VALASEK P, ROUTH P, et al. An efficient multiscale method for time-domain waveform tomography[J]. Geophysics, 2009, 74(6): WCC59-WCC68. DOI:10.1190/1.3151869 |