2. 中国石油大学(华东)地球科学与技术学院, 山东青岛 266580
2. School of Geosciences in China University of Petroleum(East China), Qingdao 266580, China
随着油气勘探开发的深入, 勘探精度要求逐渐提高, 地震波成像也逐步从构造成像向岩性成像发展[1]。岩性解释要求地震成像结果具有较高的分辨率和振幅保真度, 但常规的偏移成像方法只能提供模糊的结构影像, 且振幅属性不可靠[2]。最小二乘逆时偏移(LSRTM)方法是在最小二乘反演的理论框架下, 通过构造互为共轭的偏移算子和反偏移算子, 将逆时偏移[3]成像结果作为最小二乘反演的第一步[4], 通过不断更新成像结果, 达到为岩性解释提供高分辨率和真振幅成像的目的[5-7]。LSRTM的研究主要集中在4个方面:计算效率[8-11]、复杂介质[12-16]、目标泛函[17-19]和正则化约束[20-25]。合适的正则化约束不仅可以有效提高反演的收敛速度, 而且可以保证反演过程稳定, 改善反演效果。通过共成像点道集的聚焦性或平滑性约束, 能够压制最小二乘偏移对不完整地震数据偏移产生的成像噪音。但共成像点道集需要额外的计算量的存储空间, 因此有些学者提出将一些去噪方法(如小波变换、曲波域变换、平面波解构滤波等)引入最小二乘偏移约束反演结果的稀疏性。Dutta等[20-21]利用Seislet变换和局部Radon变换对LSRTM进行约束, 有效压制了多震源的串扰噪声, 加快了收敛速度。Dai等[22]将正则化项引入最小二乘偏移中, 加快了收敛速度, 提高了计算效率。刘玉金等[23]研究了局部倾角约束的最小二乘偏移方法, 比较了多种约束条件。王彦飞[24]分析了正则化对LSM的作用, 并介绍了干涉偏移的意义。Wu等[25]提出了L1稀疏约束的LSRTM, 有效压制了偏移噪声。当前的LSRTM大都是从常密度二阶声波方程推导的[26], 考虑到多参数反演的优势, 李庆洋等[27]推导了基于一阶速度-应力方程的相位编码LSRTM算法, 但对串扰噪音的压制不够明显, 且仅讨论了二维情况。笔者将基于一阶速度-应力方程的相位编码LSRTM算法推广到3D情形, 且为压制偏移噪音引入一种能够快速求解的稀疏正则化约束。
1 方法原理 1.1 线性Born正演方程三维声波介质下的一阶速度-应力波动方程为
$ \left\{ {\begin{array}{*{20}{l}} {\rho \frac{{\partial u}}{{\partial t}} = \frac{{\partial p}}{{\partial x}},}\\ {\rho \frac{{\partial v}}{{\partial t}} = \frac{{\partial p}}{{\partial y}},}\\ {\rho \frac{{\partial w}}{{\partial t}} = \frac{{\partial p}}{{\partial z}},}\\ {{s^2}\frac{{\partial p}}{{\partial t}} = \rho \left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}} + \frac{{\partial w}}{{\partial z}}} \right) + f.} \end{array}} \right. $ | (1) |
式中, p为声压; u、v、w分别为质点x、y、z方向速度; ρ为介质密度; s为速度的倒数, 即慢度; f为震源。
为了便于表示, 令矩阵
$ \mathit{\boldsymbol{A}} = \left( {\begin{array}{*{20}{c}} 0&0&0&{ - 1}\\ 0&0&0&0\\ 0&0&0&0\\ { - \rho }&0&0&0 \end{array}} \right),\mathit{\boldsymbol{B}} = \left( {\begin{array}{*{20}{c}} 0&0&0&0\\ 0&0&0&{ - 1}\\ 0&0&0&0\\ 0&{ - \rho }&0&0 \end{array}} \right), $ |
$ \mathit{\boldsymbol{C}} = \left( {\begin{array}{*{20}{c}} 0&0&0&0\\ 0&0&0&0\\ 0&0&0&{ - 1}\\ 0&0&{ - \rho }&0 \end{array}} \right),\mathit{\boldsymbol{D}} = \left( {\begin{array}{*{20}{c}} \rho &0&0&0\\ 0&\rho &0&0\\ 0&0&\rho &0\\ 0&0&0&{{s^2}} \end{array}} \right), $ |
$ 向量\;\mathit{\boldsymbol{U}} = \left( {\begin{array}{*{20}{l}} u\\ v\\ w\\ p \end{array}} \right), $ |
从而可将公式(1)简化为
$ \mathit{\boldsymbol{LU}} = \mathit{\boldsymbol{F}}. $ | (2) |
其中
$ \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{A}}{\partial _x} + \mathit{\boldsymbol{B}}{\partial _y} + \mathit{\boldsymbol{C}}{\partial _z} + \mathit{\boldsymbol{D}}{\partial _t}, $ |
$ \mathit{\boldsymbol{F}} = \left( {\begin{array}{*{20}{l}} 0&0&0&{f{)^{\rm{T}}}} \end{array}} \right.. $ |
式中, L为正演算子; F为震源矩阵。
LSRTM的第一步是构建线性Born正演方程。假定介质密度不变, 将总慢度场s2写成背景慢度场s02和扰动慢度场Δs2的加和形式, 即
$ {s^2} = s_0^2 + \Delta {s^2}. $ | (3) |
由场的叠加原理可知, 原始波场U可理解为由背景介质产生的背景波场U0和由扰动介质产生的扰动波场Us叠加而成, 即
$ \mathit{\boldsymbol{U}} = {\mathit{\boldsymbol{U}}_0} + {\mathit{\boldsymbol{U}}_{\rm{s}}}. $ | (4) |
原始波场和背景波场都符合波动方程, 从而将式(3)与(4)代入式(1), 并应用Born近似即可得到扰动波场Us的控制方程
$ \left\{ {\begin{array}{*{20}{l}} {\rho \frac{{\partial {u_{\rm{s}}}}}{{\partial t}} = \frac{{\partial {p_{\rm{s}}}}}{{\partial x}},}\\ {\rho \frac{{\partial {v_{\rm{s}}}}}{{\partial t}} = \frac{{\partial {p_{\rm{s}}}}}{{\partial y}},}\\ {\rho \frac{{\partial {w_{\rm{s}}}}}{{\partial t}} = \frac{{\partial {p_{\rm{s}}}}}{{\partial z}},}\\ {{s^2}\frac{{\partial {p_{\rm{s}}}}}{{\partial t}} = \rho \left( {\frac{{\partial {u_{\rm{s}}}}}{{\partial x}} + \frac{{\partial {v_{\rm{s}}}}}{{\partial y}} + \frac{{\partial {w_{\rm{s}}}}}{{\partial z}}} \right) - \Delta {s^2}\frac{{\partial {p_0}}}{{\partial t}}.} \end{array}} \right. $ | (5) |
将其表达成矩阵的形式为
$ \mathit{\boldsymbol{L}}{\mathit{\boldsymbol{U}}_{\rm{s}}} = \mathit{\boldsymbol{\widehat F}}. $ | (6) |
其中
$ \mathit{\boldsymbol{\hat F}} = {\left( {\begin{array}{*{20}{c}} 0&0&0&{ - m\frac{{\partial {p_0}}}{{\partial t}}} \end{array}} \right)^{\rm{T}}}, $ |
$ m = \Delta {s^2}, $ |
$ {\mathit{\boldsymbol{U}}_{\rm{s}}} = {\left( {\begin{array}{*{20}{l}} {{u_{\rm{s}}}}&{{v_{\rm{s}}}}&{{w_{\rm{s}}}}&{{p_{\rm{s}}}} \end{array}} \right)^{\rm{T}}}. $ |
LSRTM是在反演框架内, 通过实际观测数据求取地下介质的反射系数。其无约束项的目标泛函为
$ J(m) = \frac{1}{2}\left\| {{\mathit{\boldsymbol{U}}_{\rm{s}}} - {\mathit{\boldsymbol{U}}_{{\rm{obs}}}}} \right\|_2^2. $ | (7) |
式中, Uobs为观测记录;
一般采用梯度导引类算法求解公式(7), 梯度可表示为
$ g(m) = \left\langle {{\mathit{\boldsymbol{U}}_{\rm{s}}} - {\mathit{\boldsymbol{U}}_{{\rm{obs}}}},\frac{{\partial {\mathit{\boldsymbol{U}}_{\rm{s}}}}}{{\partial m}}} \right\rangle . $ | (8) |
其中, 〈, 〉为标量乘;
$ \mathit{\boldsymbol{L}}\frac{{\partial {\mathit{\boldsymbol{U}}_{\rm{s}}}}}{{\partial m}} = \mathit{\boldsymbol{G}}. $ | (9) |
其中
$ \mathit{\boldsymbol{G}} = {\left( {\begin{array}{*{20}{c}} 0&0&0&{ - \frac{{\partial {p_0}}}{{\partial t}}} \end{array}} \right)^{\rm{T}}}. $ |
借助式(9)和伴随状态法, 可将式(8)改写为
$ \begin{array}{l} g(m) = = \left\langle {\mathit{\boldsymbol{U}}, - {\mathit{\boldsymbol{U}}_{{\rm{obs}}}},{\mathit{\boldsymbol{L}}^{ - 1}}(\mathit{\boldsymbol{G}})} \right\rangle = \left\langle {{{\left( {{\mathit{\boldsymbol{L}}^{ - 1}}} \right)}^ * }\left( {{\mathit{\boldsymbol{U}}_{\rm{s}}} - {\mathit{\boldsymbol{U}}_{{\rm{obs}}}}} \right),} \right.\\ \left. \mathit{\boldsymbol{G}} \right\rangle . \end{array} $ | (10) |
可以看出, 若定义伴随波场U*满足方程
$ {\mathit{\boldsymbol{L}}^*}{\mathit{\boldsymbol{U}}^*} = {\mathit{\boldsymbol{U}}_{\rm{s}}} - {\mathit{\boldsymbol{U}}_{{\rm{obs}}}}. $ | (11) |
则公式(10)中标量乘的左项可直接用伴随波场U*表示, 从而梯度项可直接计算得到。式中, L*为L的伴随算子, 借助伴随状态法进行标量乘的积分展开, 然后应用分步积分、初始条件和边界条件, 可推出L*的表达式为
$ {\mathit{\boldsymbol{L}}^*} = - {\mathit{\boldsymbol{A}}^{\rm{T}}}{\partial _x} - {\mathit{\boldsymbol{B}}^{\rm{T}}}{\partial _y} - {\mathit{\boldsymbol{C}}^{\rm{T}}}{\partial _z} - {\mathit{\boldsymbol{D}}^{\rm{T}}}{\partial _t}. $ | (12) |
将伴随算子L*的表达式代入公式(11)即可得伴随方程:
$ \left\{ {\begin{array}{*{20}{l}} {\rho \frac{{\partial {p^*}}}{{\partial x}} - \rho \frac{{\partial {u^*}}}{{\partial t}} = {u_{\rm{s}}} - {u_{{\rm{obs}}}},}\\ {\rho \frac{{\partial {p^*}}}{{\partial y}} - \rho \frac{{\partial {v^*}}}{{\partial t}} = {v_{\rm{s}}} - {v_{{\rm{obs}}}},}\\ {\rho \frac{{\partial {p^*}}}{{\partial z}} - \rho \frac{{\partial {w^*}}}{{\partial t}} = {w_{\rm{s}}} - {w_{{\rm{obs}}}},}\\ {\frac{{\partial {u^*}}}{{\partial x}} + \frac{{\partial {v^*}}}{{\partial y}} + \frac{{\partial {w^*}}}{{\partial z}} - {s^2}\frac{{\partial {p^*}}}{{\partial t}} = {p_{\rm{s}}} - {p_{{\rm{obs}}}}.} \end{array}} \right. $ | (13) |
可以看出, 该伴随方程(式13)与原始波动方程(式1)相比在形式上有很大的改变, 这是一阶速度-应力方程[27]相比于二阶常密度标量方程的根本不同。
通过式(13)计算得到伴随波场后代入式(10), 则可得到梯度公式的最终表达形式:
$ g(m) = \left\langle {{\mathit{\boldsymbol{U}}^*},\mathit{\boldsymbol{G}}} \right\rangle = - {p^*}\frac{{\partial {p_0}}}{{\partial t}}. $ | (14) |
常规LSRTM算法利用最速下降法或共轭梯度法可不断迭代更新模型。
1.3 稀疏正则化约束LSRTM地震反演问题大多为不适定问题, 一般需要加入约束项, 例如常用的基于模型L2范数约束的吉洪诺夫正则化, 它假定模型参数是光滑和连续的, 对反演结果具有平滑滤波效应, 但对串扰噪音的压制不够明显。L1范数正则化可保持地下介质反射系数的稀疏性, 提高成像分辨率, 且对压制相位编码串扰噪音有较好的效果。
L1范数正则化约束的目标泛函为
$ J(m) = \frac{1}{2}\left\| {{\mathit{\boldsymbol{U}}_{\rm{s}}} - {\mathit{\boldsymbol{U}}_{{\rm{obs}}}}} \right\|_2^2 + \lambda {\left\| m \right\|_1}. $ | (15) |
考虑到L1范数是非凸的, 公式中正则化项的梯度无法直接计算, 因而梯度导引类算法不再适用。关于此类方程的常见解法主要有贪婪算法、谱梯度投影法、迭代阈值收缩法等。在大型计算时, 后两种算法更为常用, 但也各自存在问题。谱梯度投影算法计算量较大, 且实现非常复杂。迭代阈值收缩法简单易实现, 不过反演结果对阈值的选取敏感。为此, 笔者采用线性Bregman算法求解上述问题的另一变换形式:
$ \left\{ {\begin{array}{*{20}{l}} {\min \lambda {{\left\| m \right\|}_1} + \frac{1}{2}\left\| m \right\|_2^2,}\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}\frac{1}{2}\left\| {{U_s} - {U_{{\rm{obs }}}}} \right\|_2^2 \le \delta } \end{array}} \right. $ | (16) |
式中, λ为稀疏(L1范数)与平滑(L2范数)的权衡因子, 其选取方案在模型试算部分给出; δ为噪音水平的估计, 当观测数据中不含噪声时可令其为零。
更新步骤为
$ \left\{ {\begin{array}{*{20}{l}} {{z_{k + 1}} = {z_k} - {\alpha _k}{g_k},}\\ {{m_{k + 1}} = {\rm{sign}} \left( {{z_{k + 1}}} \right)\max \left( {0,\left| {{z_{k + 1}}} \right| - \lambda } \right).} \end{array}} \right. $ | (17) |
式中, zk为第k次迭代计算的中间变量, 初始值为0, 正是由于该变量的引入, 不需对模型参数m直接使用硬阈值法, 从而增加了算法的稳健性; gk为第k次迭代的L2泛函项的梯度; αk为对应的计算步长。
综上所述, 本方法的具体计算流程为:首先根据式(6)计算扰动波场, 得到与观测记录之间的残差记录; 其次利用式(13)将残差反传得到伴随波场; 然后由式(14)计算得到梯度; 最后根据式(17)不断迭代更新模型, 得到最终的反演结果。
2 模型试算 2.1 三维洼陷模型首先以三维洼陷模型为例来验证算法的正确性及有效性。速度场如图 1(a)所示, 图 1(b)为相应的反射率模型(慢度平方的扰动)。三维网格数分别为301、201和121, 3个方向的网格间距都是10 m。由于三维最小二乘逆时偏移的计算量过于庞大, 普通的个人计算机已无法满足多炮计算的需求, 因而在此利用多震源相位编码技术将所有炮组合成一个超级炮, 从而使计算需求降低到单炮的水平, 基本满足了单机计算的需要。计算参数为:每个X方向共30炮以100 m的间隔均匀分布在地表, Y方向有10炮以200 m的间隔均匀排列, 总共300炮; 共有10 201个检波器以X方向30 m和Y方向20 m的间隔均匀分布。震源为25 Hz的Ricker子波, 数值模拟的时间采样间隔为1 ms, 最大记录2.0 s。所用计算机为Dell工作站, 有64 GB内存, 32个线程, 计算主频为3.1 GHz。采用共享存储并行编程(OpenMP)技术后单次更新的计算时间约为1 h, 总共更新了150次, 共耗时接近1周。
图 2(a)、(c)分别为常规三维相位编码LSRTM算法迭代更新60次和150次时的成像结果, 随着迭代次数的不断增大, 成像结果中的低频和高频噪音都得到了一定程度的压制减弱, 且震源效应也有较好的消除, 整个剖面的振幅更加均衡, 不过最终的成像结果中仍然存在一些高频的偏移串扰噪音。图 2(b)、(d)为本文中L1范数正则化相位编码LSRTM算法第60次和150次的迭代结果。可以看出, 常规LSRTM算法即使随着迭代次数的不断增大, 也无法完全压制由相位编码引入的高频串扰噪音; 而L1范数正则化算法由于加入了稀疏约束, 反演结果可有效地压制成像结果中的低频和高频串扰噪音, 显著提升成像分辨率, 较大程度地改善成像质量。
相应的两种算法的模型残差曲线如图 3所示, 其中虚线为传统算法, 实线为本文算法。从图 3可以清楚地看出, 本文算法收敛更快, 在第100次迭代时的成像效果已经优于传统算法第160次的成像效果, 即采用本文算法所需迭代次数更少, 从而可显著地节约计算量、提高计算效率。
由于前几次迭代的成像结果完全被低频噪音和相位编码串扰覆盖, 因而从第一次开始就加入稀疏约束对最终的结果改善有限, 多次试验发现, 从第30到60次迭代开始加入稀疏约束较好。具体来说, 本文中从第51次开始加入稀疏约束, 即前50次迭代中正则化参数λ设为0, 第51次的正则化参数λ取为第50次反演结果的最大绝对值的10%, 之后保持正则化参数不变。
2.2 三维盐丘模型最后以3D SEG/EAGE Salt模型测试本文算法对复杂模型的适用性, 速度场如图 4(a)所示, 将该速度场平滑后计算得到的慢度平方的扰动场模型如图 4(b)所示。离散后的模型X和Y方向的网格数目都是338, Z方向则是201个网格, 空间网格步长都是20 m。计算参数为:每个X方向共30炮以200 m的间隔均匀分布在地表, Y方向有20炮以300 m的间隔均匀排列, 总共600炮; 共有10 201个检波器以X方向30 m和Y方向20 m的间隔均匀分布。正演模拟所用震源为主频12 Hz的雷克子波, 时间采样间隔为1.6 ms, 最大记录时间6.4 s。
同样为了减小计算量、提高计算效率, 利用相位编码技术将600炮组合为一个超级炮。图 5(a)和(b)分别为本文L1范数正则化LSRTM算法在第100次和160次的迭代结果。可以看出, 即使对复杂模型本文算法仍然能得到较好的效果, 随着迭代次数的增大, 高频的串扰噪音得到了一定的压制消除, 较大地改善了成像质量。当然由于模型过于复杂, 且迭代次数不够, 第160次的成像结果中仍然存在一些噪音, 这可以通过增加超级炮的数目进一步提升成像效果。
(1) 基于一阶速度-应力方程的LSRTM算法易于处理强非均匀变密度介质, 本文中将其推广到三维情形, 拓展了算法的适用范围, 从而更加符合实际情况。
(2) 多震源相位编码算法可有效降低计算量、提高计算效率, 在此基础上加入L1范数稀疏约束可较好地压制高频串扰噪声, 高精度地恢复地下反射系数, 显著改善最终成像质量。
(3) 线性Bergman解法引入辅助变量, 使反演解法更加稳健, 降低了方法对正则化参数的依赖, 更加适用于实际资料的处理, 且本文给出的参数选取准则具有一定的普适性, 可适用于全波形反演、层析反演等其他情形。
[1] |
李振春. 地震偏移成像技术研究现状与发展趋势[J]. 石油地球物理勘探, 2014, 49(1): 1-21. LI Zhenchun. Research status and development trends for seismic migration technology[J]. Oil Geophysical Prospecting, 2014, 49(1): 1-21. |
[2] |
NEMETH T, WU C J, SCHUSTER G T. Least-squares migration of incomplete reflection data[J]. Geophysics, 1999, 64(1): 208-221. DOI:10.1190/1.1444517 |
[3] |
李振春, 雍鹏, 黄建平, 等. 基于矢量波场分离弹性波逆时偏移成像[J]. 中国石油大学学报(自然科学版), 2016, 40(1): 42-48. LI Zhenchun, YONG Peng, HUANG Jianping, et al. Elastic wave reverse time migration based on vector wavefield separation[J]. Journal of China University of Petroleum (Edition of Natural Science), 2016, 40(1): 42-48. |
[4] |
TARANTOLA A. Linearized inversion of seismic reflection data[J]. Geophysical Prospecting, 1984, 32(6): 998-1015. DOI:10.1111/j.1365-2478.1984.tb00751.x |
[5] |
DAI W, FOWLER P, SCHUSTER G T. Multi-source least-squares reverse time migration[J]. Geophysical Prospecting, 2012, 60(4): 681-695. DOI:10.1111/j.1365-2478.2012.01092.x |
[6] |
郭振波, 李振春. 最小平方逆时偏移真振幅成像[J]. 石油地球物理勘探, 2014, 49(1): 113-120. GUO Zhenbo, LI Zhenchun. True-amplitude imaging based on least-squares reverse time migration[J]. Oil Geophysical Prospecting, 2014, 49(1): 113-120. |
[7] |
郭书娟, 马方正, 段心标, 等. 最小二乘逆时偏移成像方法的实现与应用研究[J]. 石油物探, 2015, 54(3): 301-308. GUO Shujuan, MA Fangzheng, DUAN Xinbiao, et al. Research of least-squares reverse-time migration imaging method and its application[J]. Geophysical Prospecting for Petroleum, 2015, 54(3): 301-308. DOI:10.3969/j.issn.1000-1441.2015.03.008 |
[8] |
ROMERO L A, GHIGLIA D C, OBER C C, et al. Phase encoding of shot records in prestack migration[J]. Geophysics, 2000, 65(2): 426-436. DOI:10.1190/1.1444737 |
[9] |
SCHUSTER G T, WANG X, HUANG Y, et al. Theory of multisource crosstalk reduction by phase-encoded statics[J]. Geophysical Journal International, 2011, 184(3): 1289-1303. DOI:10.1111/j.1365-246X.2010.04906.x |
[10] |
DAI W, SCHUSTER G T. Plane-wave least-squares reverse-time migration[J]. Geophysics, 2013, 78(4): S165-S177. DOI:10.1190/geo2012-0377.1 |
[11] |
李庆洋, 黄建平, 李振春, 等. 优化的多震源最小二乘逆时偏移[J]. 石油地球物理勘探, 2016, 51(2): 334-341. LI Qingyang, HUANG Jianping, LI Zhenchun, et al. Optimized multi-source least-squares reverse time migration[J]. Oil Geophysical Prospecting, 2016, 51(2): 334-341. |
[12] |
黄建平, 曹晓莉, 李振春, 等. 最小二乘逆时偏移在近地表高精度成像中的应用[J]. 石油地球物理勘探, 2014, 49(1): 107-112. HUANG Jianping, CAO Xiaoli, LI Zhenchun, et al. Least square reverse time migration in high resolution imaging of near surface[J]. Oil Geophysical Prospecting, 2014, 49(1): 107-112. |
[13] |
DUTTA G, SCHUSTER G T. Attenuation compensation for least-squares reverse time migration using the viscoacoustic-wave equation[J]. Geophysics, 2014, 79(6): S251-S262. DOI:10.1190/geo2013-0414.1 |
[14] |
李振春, 郭振波, 田坤. 黏声介质最小平方逆时偏移[J]. 地球物理学报, 2014, 57(1): 214-228. LI Zhenchun, GUO Zhenbo, TIAN Kun. Least-squares reverse time migration in visco-acoustic medium[J]. Chinese Journal of Geophysics, 2014, 57(1): 214-228. |
[15] |
ZHANG D, SCHUSTER G T. Least-squares reverse time migration of multiples[J]. Geophysics, 2013, 79(1): S11-S21. |
[16] |
刘学建, 刘伊克. 表面多次波最小二乘逆时偏移成像[J]. 地球物理学报, 2016, 59(9): 3354-3365. LIU Xuejian, LIU Yike. Least-squares reverse-time migration of surface-related multiples[J]. Chinese Journal of Geophysics, 2016, 59(9): 3354-3365. |
[17] |
TAN S, HUANG L. Least-squares reverse-time migration with a wavefield-separation imaging condition and updated source wavefields[J]. Geophysics, 2014, 79(5): S195-S205. DOI:10.1190/geo2014-0020.1 |
[18] |
LIU Y, TENG J, XU T, et al. An efficient step-length formula for correlative least-squares reverse time migration[J]. Geophysics, 2016, 81(4): S221-S238. DOI:10.1190/geo2015-0529.1 |
[19] |
ZHANG Q, ZHOU H, CHEN H, et al. Least-squares reverse time migration with and without source wavelet estimation[J]. Journal of Applied Geophysics, 2016, 134: 1-10. DOI:10.1016/j.jappgeo.2016.08.003 |
[20] |
DUTTA G. Sparse least-squares reverse time migration using seislets[J]. Journal of Applied Geophysics, 2017, 136: 142-155. DOI:10.1016/j.jappgeo.2016.10.027 |
[21] |
DUTTA G, GIBOLI M, AGUT C, et al. Least-squares reverse time migration with local Radon-based preconditioning[J]. Geophysics, 2017, 82(2): S75-S84. DOI:10.1190/geo2016-0117.1 |
[22] |
DAI W, WANG X, SCHUSTER G T. Least-squares migration of multisource data with a deblurring filter[J]. Geophysics, 2011, 76(5): R135-R146. DOI:10.1190/geo2010-0159.1 |
[23] |
刘玉金, 李振春, 吴丹, 等. 局部倾角约束最小二乘偏移方法研究[J]. 地球物理学报, 2013, 56(3): 1003-1011. LIU Yujin, LI Zhenchun, WU Dan, et al. The research on local slope constrained least-squares migration[J]. Chinese Journal of Geophysics, 2013, 56(3): 1003-1011. |
[24] |
王彦飞. 地震波干涉偏移及预条件正则化最小二乘偏移成像方法对比[J]. 地球物理学报, 2013, 56(1): 230-238. WANG Yanfei. Comparison of interferometric migration and preconditioned regularizing least squares migration inversion methods in seismic imaging[J]. Chinese Journal of Geophysics, 2013, 56(1): 230-238. |
[25] |
WU D, YAO G, CAO J, et al. Least-squares RTM with L1 norm regularization[J]. Journal of Geophysics and Engineering, 2016, 13(5): 666-673. DOI:10.1088/1742-2132/13/5/666 |
[26] |
李庆洋, 黄建平, 李振春, 等. 基于一阶速度-应力方程的多震源最小二乘逆时偏移[J]. 地球物理学报, 2016, 59(12): 4666-4676. LI Qingyang, HUANG Jianping, LI Zhenchun, et al. Multi-source least-squares reverse time migration based on first-order velocity-stress wave equation[J]. Chinese Journal of Geophysics, 2016, 59(12): 4666-4676. DOI:10.6038/cjg20161226 |
[27] |
黄建平, 杨宇, 李振春, 等. 基于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. |