• ISSN 2096-8957
  • CN 10-1702/P

地震逆散射偏移与反演综述

毛伟建 李武群 欧阳威

引用本文:
Citation:

地震逆散射偏移与反演综述

Review of seismic inverse scattering migration and inversion

    Corresponding author: Mao Weijian, wjmao@whigg.ac.cn ;
  • CLC number: P315.61, P315.63

  • 摘要: 随着油气勘探领域逐渐向深层、复杂型、隐蔽性油气藏转移,油气资源的勘探难度越来越大,传统反射地震勘探技术难以满足日益增长的油气勘探需求,亟需发展适合复杂地质构造的地震波偏移反演新技术. 针对地球深部非均匀结构体引起的地震散射波,发展地震逆散射偏移反演理论和技术将有可能解决复杂构造成像反演的技术难题. 本文回顾地震波逆散射偏移反演理论的发展历史和基本原理,以逆广义Radon变换求解线性化逆散射问题为基础,介绍逆散射理论在介质结构成像、物性参数反演、多次波衰减等方面的技术延伸,同时将其应用到合成数据和实际数据资料,探讨地震勘探逆散射方法的技术优势和应用潜力.
  • 图 1  (a)常规GRT反演剖面;(b)角度域GRT反演剖面(修改自Li et al., 2018

    Figure 1.  Inversion profiles for the layered model by using (a) conventional GRT inversion method and (b) angle-domain GRT inversion method (modified from Li et al., 2018)

    图 2  常规GRT反演振幅曲线(黑色虚线)和真实扰动振幅(红色实线)对比图. 道集位置对应图1a中的白色虚线位置: (a)x=500 m;(b)x=800 m;以及(c)x=2000 m(修改自Li et al., 2018

    Figure 2.  Comparison of true perturbations (red solid lines) and retrieved perturbations (black dashed lines) by conventional GRT inversion method at the distances (a) x=500 m ; (b) x=800 m ; and (c) x=2000 m in Fig. 1a (modified from Li et al., 2018)

    图 3  图2类似,角度域GRT反演振幅曲线与真实扰动振幅对比图(修改自Li et al., 2018

    Figure 3.  Similar to Fig. 2, but for angle-domain GRT inversion method from Fig. 1b (modified from Li et al., 2018)

    图 4  水平保幅情况.(a)常规GRT反演;(b)角度域GRT反演(修改自Li et al., 2018

    Figure 4.  The retrieved amplitude curves picked along the reconstructed reflectors from (a) conventional GRT inversion profile and (b) angle-domain GRT inversion profile (modified from Li et al., 2018)

    图 5  变密度水平层状模型角度域共成像点道集图(修改自Li et al., 2018

    Figure 5.  ADCIGs calculated at the distance x=2 000 m (modified from Li et al., 2018)

    图 6  (a)理论AVA特征曲线;(b)计算提取的AVA曲线(修改自Li et al., 2018

    Figure 6.  (a) Theoretical AVA curves and (b) Picked AVA curves from the calculated ADCIGs panel (Fig. 5) (modified from Li et al., 2018)

    图 7  散射序列示意图.(a)全波场散射示意图.(b)一阶Born近似散射示意图.(c)二阶Born近似散射示意图(李武群等, 2017

    Figure 7.  Schematics of scattering series for (a) full order Born scattering, (b) first order scattering, and (c) second order scattering (Li et al., 2017)

    图 8  局部二次散射示意图(李武群等,2017

    Figure 8.  Schematics of second order scattering within a local area (Li et al., 2017)

    图 9  不同速度扰动的单界面水平层状模型

    Figure 9.  Single interface horizontal layered model with different velocity perturbations

    图 10  不同扰动模型Born近似误差变化曲线.(a)一阶Born近似误差.(b)二阶Born近似误差(李武群等,2017

    Figure 10.  Error distribution curves from Born data of different perturbation. (a) First order Born approximate data error. (b) Second order Born approximate data error (Li et al., 2017)

    图 11  GRT线性反演和非线性反演界面振幅归一化曲线图.(a)GRT线性反演;(b)GRT非线性反演(李武群等,2017

    Figure 11.  Comparison of normalized inverted amplitudes from GRT linear and nonlinear inversion of different perturbation models. (a) GRT linear inversion. (b) GRT nonlinear inversion (Li et al., 2017)

    图 12  (a)Sigsbee 2A速度模型及(b)非线性反演剖面(修改自Ouyang et al., 2014

    Figure 12.  (a) the Sigsbee 2A model and (b) the nonlinear inversion profile (modified from Ouyang et al., 2014)

    图 13  反演值(实线)与真实扰动(虚线)之间的比较.(a)线性反演值与真实扰动之间的比较.(b)非线性反演值与真实扰动的比较(修改自Ouyang et al., 2014

    Figure 13.  Comparison between true perturbations (dotted line) and inverted perturbations (solid line). (a) Linear inversion. (b) Quadric nonlinear inversion (modified from Ouyang et al., 2014)

    图 14  墨西哥湾某区ISS自由表面多次波去除应用示例(修改自Weglein et al., 2003

    Figure 14.  The left panel is a stack of a field data set from the Gulf of Mexico. The right panel is the result of inverse-scattering free-surface multiple removal (modified from Weglein et al., 2003)

    图 15  墨西哥湾某区ISS层间多次波去除应用示例(修改自Weglein et al., 2003

    Figure 15.  An example of inverse-scattering internal multiple attenuation from the Gulf of Mexico (modified from Weglein et al., 2003)

    图 16  二维SEAM模型

    Figure 16.  The 2D SEAM model

    图 17  2D SEAM 模型偏移成像结果.(a)逆散射成像;(b)高斯束成像

    Figure 17.  The migration results of the 2D SEAM model by (a) ISM, and (b) GBM

    图 18  3D SEG/EAGE盐丘模型

    Figure 18.  The 3D SEG/EAGE salt model

    图 19  3D SEG/EAGE模型沿x方向剖面偏移成像结果.(a)和(b)分别原始速度模型.(c)和(d)为逆散射偏移成像结果.(e)和(f)为高斯束偏移成像结果

    Figure 19.  Migration results of the profiles along x direction. (a) and (b) are velocity models. (c) and (d) are ISM results. (e) and (f) are GBM results

    图 20  3D SEG/EAGE模型沿y方向剖面偏移成像结果.(a)和(b)分别原始速度模型.(c)和(d)为逆散射偏移成像结果.(e)和(f)为高斯束偏移成像结果

    Figure 20.  Migration results of the profiles along y direction. (a) and (b) are velocity models. (c) and (d) are ISM results. (e) and (f) are GBM results

    图 21  3D SEG/EAGE模型深度切片偏移成像结果. (a)为深度位置z=0.8 km切片上的速度模型. (b)为逆散射偏移成像结果. (c)为高斯束偏移成像结果

    Figure 21.  Migration results of the depth slice. (a) Velocity at slice z=0.8 km. (b) ISM result. (c) GBM result

    图 22  长排列采集观测系统

    Figure 22.  The Schematic of the acquisition geometry

    图 23  野外单炮地震记录

    Figure 23.  The single shot record of the field data

    图 24  逆散射偏移成像剖面(左)和局部放大区域(右)

    Figure 24.  The imaging result by Inverse scattering migration (left) and the zoomed one (right)

    图 25  Kirchhoff偏移成像剖面(左)和局部放大区域(右)

    Figure 25.  The imaging result by Kirchhoff migration (left) and the zoomed one (right)

    图 26  SEAM成像结果. (a),(b),(c)和(d)依次为速度场,上行波成像,下行波成像以及全波成像结果(修改自Li et al., 2019

    Figure 26.  Imaging results of the SEAM 2D model (a). (b), (c) and (d) are the imaging results using primary, multiple and joint imaging condition, respectively (modified from Li et al., 2019)

    表 1  常密度水平层状模型

    Table 1.  Horizontal layered model with constant density

    层厚度/m速度/(m·s–1)密度/(kg·m–3)
    30020002500
    30021002500
    30020002500
    30021002500
    30020002500
    30021002500
    下载: 导出CSV

    表 2  变密度水平层状模型

    Table 2.  Horizontal layered model with varying density

    层厚度/m速度/(m·s–1)密度/(kg·m–3)
    30020002000
    30021002100
    30020002000
    30021002100
    30020002000
    30021002100
    下载: 导出CSV
  • [1]

    Beylkin G. Imaging of Discontinuities in the Inverse Scattering Problem by Inversion of a Causal Generalized Radon-Transform[J]. Journal of Mathematical Physics, 1985, 26(1): 99-108. doi: 10.1063/1.526755
    [2]

    Beylkin G. The inversion problem and applications of the generalized Radon transform[J]. Communications on pure and applied mathematics, 1984, 37(5): 579-599. doi: 10.1002/cpa.3160370503
    [3]

    Beylkin G, Burridge R. Linearized inverse scattering problems in acoustics and elasticity[J]. Wave Motion, 1990, 12(1): 15-52. doi: 10.1016/0165-2125(90)90017-X
    [4]

    Bleistein N. On the imaging of reflectors in the earth[J]. Geophysics, 1987, 52(7): 931-942. doi: 10.1190/1.1442363
    [5]

    Burridge R, De Hoop M V, Miller D E, et al., Multiparameter inversion in anisotropic elastic media[J]. Geophysical Journal International, 1998, 134(3): 757-777.
    [6] 陈小宏, 刘华锋. 预测多次波的逆散射级数方法与SRME方法及比较[J]. 地球物理学进展, 2012, 27(3): 1040-1050.

    Chen X H, Liu H F. Comparison between inverse scattering series method and SRME method in free surface related multiple prediction[J]. Progress in Geophys, 2012, 27(3): 1040-1050 (in Chinese).
    [7]

    Clayton R W, Stolt R H. A Born-WKBJ inversion method for acoustic reflection data[J]. Geophysics, 1981, 46(11): 1559-1567. doi: 10.1190/1.1441162
    [8]

    Cohen J K. Bleistein N, An inverse method for determining small variations in propagation speed[J]. Siam Journal on Applied Mathematics, 1977, 32(4): 784-799. doi: 10.1137/0132066
    [9]

    Cohen J K, Bleistein N. Velocity inversion procedure for acoustic waves[J]. Geophysics, 1979, 44(6): 1077-1087. doi: 10.1190/1.1440996
    [10]

    De Hoop M V, Bleistein N, Generalized radon transform inversions for reflectivity in anisotropic elastic media[J]. Inverse Problems, 1997, 13(3): 669-690. doi: 10.1088/0266-5611/13/3/009
    [11]

    De Hoop M V, Brandsbergdahl S. Maslov asymptotic extension of generalized Radon transform inversion in anisotropic elastic media: a least-squares approach[J]. Inverse Problems, 2000, 16(3): 519-562. doi: 10.1088/0266-5611/16/3/301
    [12]

    De Hoop M V, Spencer C, Burridge R. The resolving power of seismic amplitude data: An anisotropic inversion/migration approach[J]. Geophysics, 1999, 64(3): 852-873. doi: 10.1190/1.1444595
    [13]

    Dillon P B, A comparison between kirchhoff and grt migration on vsp data[J]. Geophysical Prospecting, 1990, 38(7): 757-777. doi: 10.1111/j.1365-2478.1990.tb01873.x
    [14] 丁科. 地震逆散射理论与深度成像[D].长沙: 中南大学, 2002.

    Ding K. Seismic inverse scattering theory and depth imaging[D]. Changsha: Central South University, 2002 (in Chinese).
    [15] 董兴朋, 滕吉文, 马学英, 等. 二阶Born近似有限频率走时敏感核[J]. 地球物理学报, 2016, 59(3): 1070-1081.

    Dong X P, Teng J W, Ma X Y, et al. Finite-frequency traveltime sensitivity kernel based on second-order Born approximation[J]. Chinese Journal of Geophysics, 59(3): 1070-1081 (in Chinese).
    [16]

    Fleury C, Vasconcelos I. Imaging condition for nonlinear scattering-based imaging: Estimate of power loss in scattering[J]. Geophysics, 2012, 77(1): S1-S18.
    [17]

    French W S. Two-dimensional and three-dimensional migration of model-experiment reflection profiles[J]. Geophysics, 1974, 39(3): 265-277. doi: 10.1190/1.1440426
    [18] 符力耘. Born序列频散方程和Born-Kirchhoff传播算子[J]. 地球物理学报, 2010, 53(2): 370-384.

    Fu L Y. Born-series dispersion equation and Born-Kirchhoff propagators[J]. Chinese Journal of Geophysics, 2010, 53(2): 370-384 (in Chinese).
    [19] 符力耘, 杨慧珠, 牟永光. 非均匀介质散射问题的体积分方程数值解法[J]. 力学学报, 1998, 30(4): 488-494.

    Fu L Y, Yang H Z, Mu Y G. Volume integral equation method for scattering problems of inhomogeneous media [J]. Acta Mechanica Sinica, 1998, 30(4): 488-494 (in Chinese).
    [20]

    Ikelle L T. Using even terms of the scattering series for deghosting and multiple attenuation of ocean‐bottom cable data[J]. Geophysics, 1999, 64(2): 579-592. doi: 10.1190/1.1444565
    [21]

    Ikelle L T, Amundsen L, Yoo S. An optimization of the inverse scattering multiple attenuation method for OBS and VC data[J]. Geophysics, 2002, 67(4): 1293-1303. doi: 10.1190/1.1500392
    [22] 金德刚, 常旭, 刘伊克. 逆散射级数法预测层间多次波的算法改进及其策略 [J]. 地球物理学报, 2008, 51(4): 1209-1217.

    Jin D G, Chang X, Liu Y K. Algorithm improvement and strategy of internal multiples prediction based on inverse scattering series method [J]. Chinese Journal of Geophysics, 2008, 51(4): 1209-1217 (in Chinese).
    [23]

    Li W, Mao W, Li X, et al. Angle-domain inverse scattering migration/inversion in isotropic media[J]. Journal of Applied Geophysics, 2018, 154: 196-209. doi: 10.1016/j.jappgeo.2018.05.006
    [24]

    Li W, Mao W, Liang Q. Joint imaging with primaries and multiples of VSP data by GRT migration[J]. Seg Technical Program Expanded Abstracts, 2019,18(10): 5325-5329.
    [25] 李武群, 毛伟建, 欧阳威, 等. 二阶Born近似及其在逆散射GRT反演中的应用[J]. 地球物理学进展, 2017, 32(2): 0645-0656.

    Li W Q, Mao W J, Ouyang W, et al. Second-order Born appoximation and its application to inverse scattering GRT inversion [J]. Progress in Geophysics, 32(2): 0645-0656 (in Chinese).
    [26] 李翔, 胡天跃. 逆散射级数法去除自由表面多次波[J]. 地球物理学报, 2009, 52(6): 1633-1640.

    Li X, Hu T Y. Surface-related multiple removal with inverse scattering series method [J]. Chinese Journal of Geophysics, 2009, 52(6): 1633-1640 (in Chinese).
    [27]

    Liang Q, Ouyang W, Mao W, et al. Scattering pattern analysis and true-amplitude generalized Radon transform migration for acoustic transversely isotropic media with a vertical axis of symmetry[J]. Geophysical Prospecting, 2020,68(4):1211-1227. doi: 10.1111/1365-2478.12921
    [28]

    Liu F, Weglein A B, Innanen K A, et al. Multi-dimensional seismic imaging using the inverse scattering series[J]. Seg Technical Program Expanded Abstracts, 2006, 25(1): 3026-3030.
    [29]

    Mao W J, Ouyang W, Li X L. Nonlinear migration inversion with second-order born approximation[C]// 75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013, London, UK, 10-13 June 2013.
    [30]

    Matson K H. The relationship between scattering theory and the primaries and multiples of reflection seismic data[J]. Journal of Seismic Exploration, 1996, 5(1): 63-78.
    [31]

    Miller D, Oristaglio M, Beylkin G. A new slant on seismic imaging: Migration and integral geometry[J]. Geophysics, 1987, 52(7): 943-964. doi: 10.1190/1.1442364
    [32]

    Moses H E. Calculation of the scattering potential from reflection coefficients[J]. Physical Review, 1956, 102(2): 559-567. doi: 10.1103/PhysRev.102.559
    [33]

    Ouyang W, Mao W. Approximate nonlinear multiparameter inversion for multicomponent single and double P-wave scattering in isotropic elastic media[J]. Geophysical Journal International, 2018, 214(1): 14-35. doi: 10.1093/gji/ggy106
    [34]

    Ouyang W, Mao W, Li W, et al. Approximate non-linear multiparameter inversion with single and double scattering seismic wavefields in acoustic media[J]. Geophysical Journal International, 2017, 208(2): 767-789. doi: 10.1093/gji/ggw422
    [35]

    Ouyang W, Mao W, Li X. Approximate solution of nonlinear double-scattering inversion for true amplitude imaging[J]. Geophysics, 2014, 80(1): R43-R55.
    [36]

    Prosser R T. Formal solutions of inverse scattering problems. III[J]. Journal of Mathematical Physics, 1969, 21(10): 2648-2653.
    [37]

    Protasov M, Tcheverda V. True/preserving amplitude seismic imaging based on Gaussian beams application[J]. Seg Technical Program Expanded Abstracts, 2006,25(1), doi:10.1190/1.2369957.
    [38]

    Protasov M I, Tcheverda V A. True amplitude elastic Gaussian beam imaging of multicomponent walkaway vertical seismic profiling data[J]. Geophysical Prospecting, 2012, 60(6): 1030-1042. doi: 10.1111/j.1365-2478.2012.01068.x
    [39]

    Protasov M L. Tcheverda V A. True amplitude imaging by inverse generalized Radon transform based on Gaussian beam decomposition of the acoustic Green's function[J]. Geophysical Prospecting, 2011, 59(2): 197-209. doi: 10.1111/j.1365-2478.2010.00920.x
    [40]

    Ravasi M. Curtis A. Nonlinear scattering based imaging in elastic media: Theory, theorems, and imaging conditions[J]. Geophysics, 2013, 78(3): S137-S155. doi: 10.1190/geo2012-0286.1
    [41]

    Raz S. Direct reconstruction of velocity and density profiles from scattered field data[J]. Geophysics, 1981, 46(6): 832-836. doi: 10.1190/1.1441220
    [42]

    Thorbecke J W. Draganov D. Finite-difference modeling experiments for seismic interferometry[J]. Geophysics, 2011, 76(6): H1-H18. doi: 10.1190/geo2010-0039.1
    [43]

    Weglein A B. How can the inverse‐scattering method really predict and subtract all multiples from a multidimensional earth with absolutely no subsurface information?[J]. Geophysics, 1999, 18(1): 132-136.
    [44]

    Weglein A B, Araujo F V, Carvalho P M, et al. Inverse scattering series and seismic exploration[J]. Inverse Problems, 2003, 19 (6):R27-R83. doi: 10.1088/0266-5611/19/6/R01
    [45]

    Weglein A B, Gasparotto F A, Carvalho P M, et al. An inverse‐scattering series method for attenuating multiples in seismic reflection data[J]. Geophysics, 1997, 62(6): 1975-1989. doi: 10.1190/1.1444298
    [46]

    Weglein A B. Gray S H. The sensitivity of Born inversion to the choice of reference velocity; a simple example[J]. Geophysics, 1983, 48(1): 36-38. doi: 10.1190/1.1441404
    [47] 杨金龙, 朱立华. 逆散射级数层间多次波压制方法及其应用[J]. 石油物探, 2018, 57(6): 853-861.

    Yang J L, Zhu L H. Inverse scattering series internal multiple attenuation method and its application[J]. Geophysical Prospecting for Petroleum, 2018, 57(6): 853-861 (in Chinese).
    [48]

    Zhang H. Weglein A B. Direct nonlinear inversion of 1D acoustic media using inverse scattering subseries[J]. Geophysics, 2009a, 74(6): WCD29-WCD39. doi: 10.1190/1.3256283
    [49]

    Zhang H, Weglein A B. Direct nonlinear inversion of multiparameter 1D elastic media using the inverse scattering series[J]. Geophysics, 2009b, 74(6): WCD15-WCD27. doi: 10.1190/1.3251271
  • [1] V. Vavryčuk高熹微宋程万永革吕春来 . 根据震源机制迭代联合反演应力和断层取向. 地球与行星物理论评, 2014, 45(4): 21-33.
    [2] Jiajun ChongSidao NiRisheng ChuPaul Somerville张尧田宝峰 . 体波接收函数与瑞利波椭圆率联合反演研究. 地球与行星物理论评, 2017, 48(3): 218-233. doi: 10.16738/j.cnki.issn.1003-3238.201703003
    [3] V. VavryukD. Kühn赵仲和吕春来 . 波形矩张量反演的时间一频率两步法. 地球与行星物理论评, 2012, 43(6): 28-48.
    [4] P. Martínez-GarzónY. Ben-ZionN. AbolfathianG. KwiatekM. Bohnhoff李振月李泽潇万永革 . 基于震源机制解反演应力场的一种改进方法. 地球与行星物理论评, 2020, 51(2): 161-184. doi: 10.16738/j.cnki.issn.1003-3238.202002004
    [5] J. L. LiH. J. ZhangH. S. KuleliM. N. Toksoz许鑫李振月万永革 . 高频波形匹配反演震源机制解及其在小震级诱发地震中的应用. 地球与行星物理论评, 2020, 51(5): 494-508.
    [6] P. S. EarleS. RostP. M. ShearerC. Thomas林秀娜曲保安吕春来 . 近距离观测的散射P'P'波. 地球与行星物理论评, 2012, 43(4): 47-61.
    [7] Yasuko OkuyamaTakahiro FunatsuTakashi FujiiNaohiko TakamotoToshiyuki Tosha仵柯田孙凤霞杜建国朱玉萍 . 与日本中部以北松代震群(1965~1967年)相关的地壳中部流体:地球化学反演. 地球与行星物理论评, 2017, 48(1): 13-31. doi: 10.16738/j.cnki.issn.1003-3238.201701002
    [8] C. SinghM. ShekarA. SinghR. K. Chadha张尧徐沁吕春来 . 根据Lg波Q值反演得到的Hi-CLIMB项目西藏剖面地震衰减特征. 地球与行星物理论评, 2012, 43(5): 12-19.
    [9] Haruo SatoMichael C. FehlerTakuto Maeda朱维吴何珍 . 《非均匀地球中的地震波传播和散射:第二版》引言. 地球与行星物理论评, 2016, 47(4): 269-281. doi: 10.16738/j.cnki.issn.1003-3238.201604001
    [10] X. LiuD. P. Zhao张晓曼赵小艳 . 利用地方震、远震走时和面波数据联合反演日本俯冲带P波和S波层析成像. 地球与行星物理论评, 2019, 50(1): 35-63. doi: 10.16738/j.cnki.issn.1003-3238.201901003
    [11] Yongge WanShuzhong ShengJichao HuangXiang LiXin Chen杨帆盛书中朱玉萍 . 基于震源机制解数据反演构造应力张量的网格搜索法及其在中国、越南和老挝边界地区的应用. 地球与行星物理论评, 2017, 48(2): 169-184. doi: 10.16738/j.cnki.issn.1003-3238.201702004
    [12] Xuewei BaoXiaoxiao SunMingjie XuDavid W. EatonXiaodong SongLiangshu WangZhifeng DingNing MiHua LiDayong YuZhouchuan HuangPan Wang孙晓晓鲍学伟徐鸣洁朱玉萍 . 通过接收函数和瑞利波联合反演揭示青藏高原东南缘两个壳内低速通道. 地球与行星物理论评, 2016, 47(4): 329-343. doi: 10.16738/j.cnki.issn.1003-3238.201604005
    [13] E. TrasattiC. KyriakopoulosM. Chini李守勇吕春来 . 用有限元反演2009年拉奎拉(意大利)MW6.3地震的DInSAR数据. 地球与行星物理论评, 2012, 43(1): 43-49.
    [14] P. Martínez-GarzónG. KwiatekM. IckrathM. Bohnhoff崔华伟万永革吕春来 . MSATSI:结合可靠经典方法的新简化用户处理及可视化工具的应力反演MATLAB软件包. 地球与行星物理论评, 2014, 45(4): 34-45.
    [15] S. HartzellC. MendozaL. Ramirez-GuzmanY. H. ZengW. Mooney李万金赵仲和 . 中国汶川2008年MW7.9地震的破裂历史:大地测量、远震和强震动数据的单独及联合反演评价. 地球与行星物理论评, 2018, 49(3): 222-242. doi: 10.16738/j.cnki.issn.1003-3238.201803003
    [16] Y. -H. YangJ. -C. HuQ. ChenZ. -G. WangM. -C. Tsai陈静杨莹辉 . 中国北天山褶皱逆冲带2016年MW6.0呼图壁地震-盲逆冲和上覆褶皱地震. 地球与行星物理论评, 2020, 51(3): 261-272. doi: 10.16738/j.cnki.issn.1003-3238.202003005
    [17] Jeanne L. Hardebeck靳志同万永革李一琼 . 俯冲带的应力方向和俯冲带大型逆断层的强度. 地球与行星物理论评, 2016, 47(3): 255-260. doi: 10.16738/j.cnki.issn.1003-3238.201603007
    [18] P. ShearerR. Bürgmann赵雯佳吕春来 . 从2004年苏门答腊-安达曼巨型逆冲断层破裂得到的经验教训. 地球与行星物理论评, 2012, 43(1): 1-23.
    [19] E. T. BakerD. A. TennantR. A. FeelyG. T. LebonS. L. Walker陈升陶春辉 . 从野外和室内实验室两方面研究热液羽状流中颗粒物尺寸和成分对光学后向散射测量的影响. 地球与行星物理论评, 2019, 50(2): 190-200. doi: 10.16738/j.cnki.issn.1003-3238.201902005
    [20] L. BaiS. L. KlempererJ. MoriM. S. KarplusL. DingH. B. LiuG. H. LiB. W. SongS. Dhakal江勇李国辉白玲 . 主喜马拉雅逆冲断裂的横向变化控制2015年尼泊尔廓尔喀地震的破裂长度. 地球与行星物理论评, 2020, 51(1): 22-32. doi: 10.16738/j.cnki.issn.1003-3238.202001002
  • 加载中
图(26) / 表(2)
计量
  • 文章访问数:  741
  • HTML全文浏览量:  352
  • PDF下载量:  23
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-07-17
  • 网络出版日期:  2020-09-09
  • 刊出日期:  2021-01-01

地震逆散射偏移与反演综述

摘要: 随着油气勘探领域逐渐向深层、复杂型、隐蔽性油气藏转移,油气资源的勘探难度越来越大,传统反射地震勘探技术难以满足日益增长的油气勘探需求,亟需发展适合复杂地质构造的地震波偏移反演新技术. 针对地球深部非均匀结构体引起的地震散射波,发展地震逆散射偏移反演理论和技术将有可能解决复杂构造成像反演的技术难题. 本文回顾地震波逆散射偏移反演理论的发展历史和基本原理,以逆广义Radon变换求解线性化逆散射问题为基础,介绍逆散射理论在介质结构成像、物性参数反演、多次波衰减等方面的技术延伸,同时将其应用到合成数据和实际数据资料,探讨地震勘探逆散射方法的技术优势和应用潜力.

English Abstract

    • 地震勘探的任务是根据记录到的波场信息确定地下介质的物性参数;通过估计岩石和流体特征的类型和范围来定位和描述地下目标体,以评估它们含油气的潜力等. 多次覆盖反射波地震勘探技术已经非常成熟,能够较好地适应横向速度变化较平缓的水平层状或小倾角构造成像问题,在浅层油气资源探测中发挥了重要的作用. 然而,面对地质结构复杂的地下深层介质,其主要表现为地层陡峭、断裂发育、小断块及破碎带等不均匀构造较多,传统的反射地震勘探技术难以满足油气勘探需求. 因此,迫切需要发展更完善的、适合复杂地质结构的地震波成像和反演新技术.

      研究表明,针对地下相对于地震波长较小的局部非均匀体,地震散射波理论更具普适性. 地震散射波的产生机理满足惠更斯—菲涅尔原理,能够从微观层面解释地震波反射、折射、绕射等单次和多次散射现象,能够更充分地利用地震波所携带的地下介质物性信息. 如何利用散射波场来重构地下介质的几何结构和物性参数是勘探地球物理界的热点问题之一. 地震逆散射成像反演理论的研究将为克服复杂构造成像技术难点提供可能. 该方法源于量子力学领域,通过研究散射体外部接收到的各种场来估计其内部结构信息. 广义而言,地震数据的处理解释、医学超声成像、无损探测等问题都可以归结为逆散射问题. 1970年代末期开始,随着波动方程逆散射研究的不断深入,该理论被引入到地球物理勘探领域并得到逐渐发展.

      经典地震逆散射理论包含两个方面:基于扰动理论的线性化散射波场正演和求解线性逆散射反问题. Cohen和Bleistein(1977, 1979)对常背景和变背景速度下的波速扰动反演进行了系统的研究,建立起了以小扰动理论和Born近似为基础的速度反演理论,并将他们所提出的方法应用到二维声波速度模型中. Raz(1981)在层状双参数声波介质模型中考虑了常数背景下速度和密度直接线性反演方法. Clayton和Stolt(1981)利用Green函数的WKBJ近似将Cohen和Bleistein、Raz的结果推广到三维变背景声波介质情形,在频率波数域中推导出速度和密度线性解. 从应用的观点出发,勘探物理学家希望知道物性参数发生突变或者间断的位置,及其梯度变化的大小(真振幅成像和反演). Miller等(1987)提出了线性真振幅成像和反演的最初轮廓−求解广义Radon变换(Generalized Radon Transform, GRT)的反投影算子,进一步完善了绕射叠加几何方法(French, 1974),而且更适于处理复杂地质构造及震源和检波器任意排列的情况. Beylkin(1984, 1985)利用小扰动技术和Born近似将反问题线性化,引入了一种Fourier积分算子并忽略其所有光滑项(只保留奇性项),有效解决了具有震荡积分核的第一类Fredholm积分方程的求解问题,奠定了奇性反演(即间断面反演)的理论基础. 在此基础上将GRT方法推广到各向同性弹性介质中,提出了与单散射夹角有关的多参数反演方法(Beylkin and Burridge, 1990). De Hoop等在各向异性介质中利用Kirchhoff近似以及逆GRT讨论了反射系数的成像和反演(De Hoop and Bleistein, 1997; De Hoop et al., 1999),而Burridge等(1998)在各向异性介质中利用Born近似和逆GRT得到弹性参数真振幅反演公式. De Hoop和Brandsbergdahl(2000)在各向异性介质中存在焦散的情况下推导了线性化的近似反演解(最小二乘的意义下),其中关键利用了Maslov近似算子来描述背景波场,对高维稳相法进行了精细分析. 从Beylkin和Miller线性逆散射GRT保幅偏移成像的观点出发,Protasov和Tcheverda(2006, 2011, 2012)推导了基于高斯束的线性逆GRT保幅偏移成像公式,特别给出了如何结合所有炮点激发的地震波信息进行角度域高斯束真振幅偏移的方法. Li等(2018)将基于GRT的线性逆散射偏移反演方法扩展到角度域,提出角度域GRT反演方法可以有效规避反演的病态问题,获得更准确的参数估计. Liang等(2020)基于Born近似和椭圆各项异性背景,提出了利用GRT近似反演声VTI介质扰动参数的方法.

      Miller最初提出的基于高频近似和Born近似的线性化逆散射GRT偏移反演理论,因其高效、灵活、保幅以及适用于任意观测系统等特性,成为非线性方程进行线性化处理的重要手段,在地震数据处理中起到了重要作用,并得到很好的发展与推广. 目前大多数偏移成像和反演算法都基于这一单散射的假设条件. 这是因为线性化的逆散射问题要比非线性逆散射问题更容易处理. 但必须承认的是,线性化的假设需要满足较为苛刻的弱散射条件,也即介质的扰动不能太大. 当地下目标区域的几何结构较为复杂,速度变化较剧烈的时候,地质条件已经不满足单散射假设,此时利用线性化逆散射偏移反演方法估计的地下介质结构物理参数并不准确.

      Born近似的适用性依赖于背景格林函数算子模拟两点之间的直达波的准确度,如果背景格林算子描述直达波不够准确,这种情况下Born近似值得怀疑. 另外,Born近似省略了序列中的高阶项,忽略了非线性多散射效应对成像的作用,这种近似处理使得线性化偏移在面对复杂介质反演振幅时难以保真. 只有将散射序列中的高阶项用于反演框架中,才能从本质上克服线性小扰动假设的不足. Moses(1956)Prosser(1969)等最早提出了逆散射级数方法(Inverse Scattering Series, ISS). 1980年代开始,研究学者为了使逆散射理论适用于大扰动的情形,开始关注逆散射序列中的高阶项. Clayton和Stolt(1981)表明通过利用散射序列中的高次项可以获得高阶近似. 勘探地球物理学家开始从散射级数角度出发研究非线性逆散射偏移反演技术. Weglein等(1999, 2003)在一维声波单参数介质中对主波场分别利用逆散射成像子级数和逆散射反演子级数进行成像和反演. Liu等(2006)将ISS成像方法推广到二维声波单参数介质中. Zhang和Weglein(2009a, 2009b)分别在一维多参数声波介质和一维各向同性弹性介质中考虑了ISS反演方法. 从ISS成像技术的最新进展可以看出该方法的数值试验都是在二维简单声波模型中进行,没有用复杂模型检测其有效性,并且假定所有多次波都已消除. 多次波去除是ISS方法一个重要的应用方向. 勘探学者1990年代中期开始研究逆散射序列,利用逆散射理论去除地震数据中的多次波(Matson, 1996; Weglein et al., 1997; Ikelle, 1999; Ikelle et al., 2002). 李翔和胡天跃(2009)将能量最低法则、带限PRT法与逆散射级数法相结合,实现了该方法在合成与实际资料中对自由表面多次波的有效压制,并详细讨论了其应用流程. 金德刚等(2008)改进了Weglein提出的ISS法预测层间多次波的1-D公式,推导了1.5-D时间—空间域ISS算法,并讨论了ISS方法在实际应用中的特点. 杨金龙和朱立华(2018)改进常规层间多次波压制处理流程,在逆散射级数法预测层间多次波前、后去除和补偿子波来提高层间多次波预测的准确性,降低对自适应相减的依赖. 陈小宏和刘华峰(2012)比较了1.5-D、2-D逆散射级数方法、和2-D SRME方法预测自由表面多次波的效果及优缺点. 除了ISS研究之外,Mao等(2013)提出了利用二阶Born近似作为一个改正项来处理单参数声波介质中的非线性多散射问题. 李武群等(2017)通过数值模拟分析了一阶Born 近似和二阶Born 近似模拟波场的准确度,并验证了二次散射局部圆域的假设. Ouyang等人对基于二阶Born 近似的非线性逆散射反演问题进行了系统深入的研究,推导了声波单参数、双参数、弹性等介质中的非线性反演公式并通过模型试算得到验证(Ouyang et al., 2014, 2017; Ouyang and Mao, 2018). 多散射问题是一个崭新的、具有挑战性的研究方向. 关于多散射的文献还可以参考丁科(2002)Fleury和Vasconcelos(2012)Ravasi和Curtis(2013)董兴朋等(2016).

      散射波理论作为地震波传播的最一般的理论,从更广义的角度去描述地震波信号,在解决复杂地质问题上具有更广泛的意义. 理论上能建立波场与地下介质扰动的数学关系,为实现高精度、高分辨率以及高保幅的地震偏移或反演提供了可能. 本文首先介绍地震波散射的基本理论和基于GRT的逆散射偏移反演方法的基本原理,围绕该方法,分别介绍角度域逆散射反演,基于二阶Born近似的非线性反演,逆散射级数预测多次波等方法. 通过合成数据和实际资料分别介绍逆散射方法在反演和成像中的典型应用和新近进展.

    • 地震波散射理论将地球介质视为背景介质和扰动介质的叠加. 震源激发的地震波在背景介质中传播形成背景波场,遇到扰动介质,背景波场与扰动介质发生一次散射形成新的震源,新的震源形成新的波前子震源继续传播,与其他扰动介质发生多次散射作用,形成更高阶散射波. 这些不同阶次的散射波最终在地表观测点上相互叠加形成非线性地震波场,它包含了地下介质的所有信息. 因此散射理论建立了介质扰动与散射波场的非线性关系. 如何通过散射波场求解介质的扰动是一个逆散射问题.下面从两个方面简要概括散射与逆散射问题:(1)正散射是给定背景介质和扰动,输出散射波场;(2)逆散射是给定背景介质以及地表观测波场,来逆向求取地下介质的扰动.

    • 从点源常密度声波波动方程出发有:

      $ {\nabla ^2}U\left({{{x}},{{s}},{\rm{\omega }}} \right) + \dfrac{{{{\rm{\omega }}^2}}}{{{c^2}\left({{x}} \right)}}U\left({{{x}},{{s}},{\rm{\omega }}} \right) = \delta \left({{{x}} - {{s}}} \right) $

      (1)

      式中,${\nabla ^2}$是Laplace微分算子,$U$是声压场,x为地下介质中散射点坐标位置矢量,s为震源点坐标位置矢量,$c$是介质中的声波波速,$\omega $是角频率,$\delta $是震源脉冲函数. 定义$f$为散射位势,表征散射点的扰动强度:

      $f\left({{x}} \right) = {c^{ - 2}}\left({{x}} \right) - {c_{0}^{ - 2}}\left({{x}} \right)$

      (2)

      式中,${c_0}$为光滑背景速度,地震波在该介质中传播产生背景散射波场,根据以上定义并类比量子力学Lippman-Schwinger方程,地震波场$U$可以表示为:

      $ U\left({{{r}},{{s}},{\rm{\omega }}} \right) \!=\! {G_0}\left({{{r}},{{s}},{\rm{\omega }}} \right) + {{\rm{\omega }}^2}\!\int\! {{\rm{d}}{{x}} {G_0}\left({{{r}},{{x}},{\rm{\omega }}} \right)f\left({{x}} \right) } U\left({{{x}},{{s}},{\rm{\omega }}} \right)\; $

      (3)

      式中,r表示检波点位置矢量,$f$为式(2)中的散射位势,${G_0}$是背景格林函数(与${c_0}$对应),满足如下方程:

      ${\nabla ^2}{G_0}\left({{{x}},{{s}},{\rm{\omega }}} \right) + \dfrac{{{{\rm{\omega }}^2}}}{{{c^2}\left({{x}} \right)}}{G_0}\left({{{x}},{{s}},{\rm{\omega }}} \right) = \delta \left({{{x}} - {{s}}} \right)$

      (4)

      根据Weglein等(2003),用L${L_0}$表示地震波在实际和背景介质中的传播算子:

      $L: = {\nabla ^2} + \dfrac{{{{\rm{\omega }}^2}}}{{{c^2}\left({{x}} \right)}},\;\;{L_0}: = {\nabla ^2} + \dfrac{{{{\rm{\omega }}^2}}}{{c_{0}^2}\left({{x}} \right)} $

      (5)

      并定义相应的扰动算子$V{\rm{ = }}L - {L_0}$,表示标量波在传播中产生的扰动. 通过式(3)不断迭代右端非线性项$U$,可以得到无穷序列:

      $U = {G_0} + {G_0}V{G_0} + {G_0}V{G_0}V{G_0} + \cdot \cdot \cdot + {G_0}{\left({V{G_0}} \right)^n} + \cdot \cdot \cdot $

      (6)

      或者:

      ${U_s} = {\left({{U_s}} \right)_1} + {\left({{U_s}} \right)_2} + {\left({{U_s}} \right)_3} + \cdot \cdot \cdot + {\left({{U_s}} \right)_n} + \cdot \cdot \cdot $

      (7)

      式中,${U_s}$表示减去背景场${G_0}$后的散射场,${\left({{U_s}} \right)_n}: = {G_0}{\left({V{G_0}} \right)^n}$表示相应的n次散射积分项($n \geqslant 1$).

      表达式(7)是散射场的完整表示形式,如果只保留其中的首阶项${\left({{U_s}} \right)_1}$,则得到散射场的一阶Born近似:

      ${U_s}\left({{{r}},{{s}},{\rm{\omega }}} \right) \approx {\left({{U_s}} \right)_1} = {{\rm{\omega }}^2}\int {{\rm{d}}{{x}} {G_0}\left({{{r}},{{x}},{\rm{\omega }}} \right)f\left({{x}} \right)} {G_0}\left({{{x}},{{s}},{\rm{\omega }}} \right)$

      (8)

      在高频近似条件下,格林函数(三维)可以表示为:

      $ {G_0}\left({{{y}},{{x}},{\rm{\omega }}} \right) \approx A\left({{{y}},{{x}}} \right)\exp \left({{{i}}{\rm{\omega }}\tau \left({{{y}},{{x}}} \right)} \right) $

      (9)

      式中,旅行时(或相位)函数$\tau \left({{{y}},{{x}}} \right)$满足程函方程:

      ${\left| {{\nabla _{{x}}}\tau \left({{{y}},{{x}}} \right)} \right|^2} = {c_0^{ - 2}}\left({{x}} \right)$

      (10)

      传播振幅$A\left({{{y}},{{x}}} \right)$满足输运方程:

      $A\left({{{y}},{{x}}} \right)\nabla _x^2\tau \left({{{y}},{{x}}} \right) + 2\nabla _x^{}A\left({{{y}},{{x}}} \right) \cdot {\nabla _x}\tau \left({{{y}},{{x}}} \right) = 0$

      (11)

      将式(9)代入式(8),并利用傅里叶变换得到时间域的散射场Born近似积分方程:

      ${U_s}\left({{{s}},{{r}},t} \right) \approx - \dfrac{{{{{{\partial}}} ^2}}}{{\partial {t}^2}}\int {{\rm{d}}{{x}} A\left({{{s}},{{x}},{{r}}} \right)\delta \left[ {t - \tau \left({{{s}},{{x}},{{r}}} \right)} \right]f\left({{x}} \right)} $

      (12)

      式中,$A\left({{{s}},{{x}},{{r}}} \right)$$\tau \left({{{s}},{{x}},{{r}}} \right)$分别为总传播振幅和旅行时函数,且有以下关系:

      $A\left({{{s}},{{x}},{{r}}} \right) = A\left({{{s}},{{x}}} \right)A\left({{{x}},{{r}}} \right), \;\; \tau \left({{{s}},{{x}},{{r}}} \right) = \tau \left({{{s}},{{x}}} \right) + \tau \left({{{x}},{{r}}} \right)$

      (13)

      散射场一阶Born近似(或单散射近似)符合传统地震反射勘探的一次波假设,在满足小尺度弱散射条件下,单散射近似与偏移之间存在理想的线性关系. 研究学者提出并论证了直接线性反演框架,在线性化的基础上解决了逆散射反演问题. 根据Beylkin和Burridge(1990),对于任意观测系统,线性逆散射问题可以类比GRT求解. 在成像点${{z}}$处,散射位势估计$\left\langle {f\left({{z}} \right)} \right\rangle $遵循以下形式:

      $\left({{{\cal{R}}}_{}^ * {U_s}} \right)\left({{z}} \right) \approx a\left({{z}} \right)\left\langle {f\left({{z}} \right)} \right\rangle .$

      (14)

      式中,${{\cal{R}}}_{}^ * $为逆GRT保幅偏移算子,具体表达式为:

      $\left({{{\cal{R}}}_{}^ * {U_s}} \right)\left({{z}} \right) = \dfrac{1}{{{\text{π} ^2}}}\iint_{} {{\rm{d}}{{s}}} {\rm{d}}{{r}} B\left({{{s}},{{z}},{{r}}} \right){U_s}\left({{{s}},{{r}},{{t}}} \right)\left| {_{t = \tau \left({{{s}},{{z}},{{r}}} \right)}^{_{}^{}}} \right.$

      (15)

      其中:

      $B\left({{{s}},{{z}},{{r}}} \right){\rm{ = }}\dfrac{{{{\cos }^2}\left({\dfrac{\theta }{2}} \right)b\left({{{z}},\theta } \right)J\left({{{s}},{{z}}} \right)J\left({{{z}},{{r}}} \right)}}{{c_0^3\left({{z}} \right)A\left({{{s}},{{z}},{{r}}} \right)}}$

      (16)

      式中,B为反演权重因子,$b$为关于成像点${{z}}$和散射夹角$\theta $的函数,$J\left({{{s}},{{z}}} \right)$$J\left({{{z}},{{r}}} \right)$分别为震源点${{s}}$和检波点${{r}}$到成像点的坐标转换雅各比参数,其计算与射线几何扩散以及入射倾角有关. $a\left({{z}} \right)$为关于成像点处散射夹角的积分:

      $a\left({{z}} \right) = \int_{{E_\theta }\left({{z}} \right)} {{\rm{d}}\theta } b\left({{{z}},\theta } \right)\sin \left({{\theta / 2}} \right)$

      (17)

      式中,${E_\theta }\left({{z}} \right)$为成像点处散射夹角变化范围. 在单参数反演中,$a\left({{z}} \right)$为成像点处一个积分值,在多参数反演框架中则是一个照明矩阵. 根据式(14),对散射场数据应用逆GRT偏移算子后,进一步计算每个成像点处的角度积分值,即可获得线性化逆散射问题的参数近似求解或重建$\left\langle {f\left({{z}} \right)} \right\rangle $.

    • 前文基于GRT的线性化逆散射问题的求解从某种意义上看,是一种“叠后”反演. 其本质是,在加权叠加后,根据成像点处的射线照明进行真振幅补偿. 因此反演的准确度非常依赖照明矩阵的计算. 该矩阵的元素是关于散射夹角的积分值,其积分限是基于成像点处的最大散射夹角(上限)和最小散射夹角(下限)来确定,而且假设积分限中夹角分布均匀. 在实际情况中,存在观测系统覆盖不充分,或者遇到陡倾结构时,成像点处的照明不充分的情形,这样确定的积分限并不准确. 在偏移过程中,只有偏移倾角与地层倾角比较接近的射线对才会对最终的偏移振幅有稳相贡献. 在照明补偿中,应该考虑的是这些有效射线对形成的夹角. 另外在多参数反演框架中,照明矩阵本身存在稳定性问题.

      Li等(2018)提出在角度域中进行反演避开常规反演的不准确、不稳定问题,发展了角度域逆散射反演方法. 该方法不将辐射模式$f\left({{{x}},\theta } \right)$中的角度信息作为照明值或矩阵独立出来,而将其作为整体进行反演. 生成角度域共成像点道集(ADCIG)后,根据成像点处具体的射线照明情况,观察叠前AVA曲线的变化特征,选择合适范围进行反演. 忽略具体过程,2D声波常密度角度域GRT偏移反演公式如下:

      $f\left({{{z}},{\theta _0}} \right) = \dfrac{1}{{\Delta \theta }}\iint_{} { {\rm{d}}{{s}}} {\rm{d}}{{r}} B\left({{{s}},{{z}},{{r}}} \right)\delta \left({\theta - {\theta _0}} \right){U_s}\left({{{s}},{{r}},t} \right)\left| {_{t = \tau \left({{{s}},{{z}},{{r}}} \right)}^{_{}^{}}} \right.$

      (18)

      式中:

      $B\left({{{s}},{{z}},{{r}}} \right) = \dfrac{{{{\cos }^2}\dfrac{\theta }{2}J\left({{{s}},{{z}}} \right)J\left({{{r}},{{z}}} \right)}}{{{{c_0^2}}\left({{z}} \right){\kappa ^0}({{z}})A\left({{{s}},{{z}},{{r}}} \right)}} $

      (19)

      在式(18)中,狄拉克$\delta $函数起到了“门”函数的作用,即将$\left[ {{\theta _0} - {{\Delta \theta } / 2},{\theta _0}{\rm{ + }}{{\Delta \theta } / 2}} \right]$内的偏移道集都归约到某一常数倾角${\theta _0}$共成像点道集中. 角度采样间隔$\Delta \theta $的选择与观测系统参数(比如炮/道间距,频率成分等)有关. 其他各符号的意义可以参考前文. 传统的逆散射GRT偏移反演方法是先加权叠加所有的射线振幅后再通过计算成像点处的照明矩阵进行振幅补偿,从而实现参数反演. 相对地,角度域逆散射保幅偏移方法则更像一种“叠前”的反演方法,首先获得高保幅并且携带了准确AVA特征信息的ADCIG,在这一基础上,根据振幅辐射模式,利用一个最小二乘拟合算法实现相应扰动参数的反演. 该方法可以避开常规逆散射GRT反演方法中不准确、不稳定的照明矩阵带来的反演问题.

      在声波介质中,线性化逆散射中的振幅辐射模式可以表示为:

      $f\left({{{x}},\theta } \right) = {f_1}\left({{x}} \right) + {f_2}\left({{x}} \right)\cos \theta $

      (20)

      式中,${f_1} = {{{K_0}} / K} - 1$表示与体模相关的扰动,${f_2} = {{{\rho _0}}/ \rho } - 1$表示与密度相关的扰动. 对于常密度单参数情形有${f_2}{\rm{ = }}0$,此时从角度域共成像点道集提取出来的AVA曲线不随角度变化,在合理照明范围内,偏移振幅曲线呈水平状态,幅值大小即为${f_1}$或者其比例. 对于变速度和变密度双参数情形,此时偏移振幅${f_{}}$表现出角度依赖,若以$\cos \theta $为横轴,成像点处的AVA振幅曲线在合理照明范围内应该接近为截距为${f_1}$,斜率为${f_2}$的直线. 通过偏移振幅辐射模式的特点,设计合适的线性拟合策略,可以比较容易地获得扰动参数${f_l}$($l = 1, 2$)的反演.

      表1是一个水平层状速度模型参数表. 该模型含有6个水平层和5个界面,水平方向和深度方向网格数为401×361,网格间距为10 m×5 m. 各层内速度均匀分布,单参数情况,各层密度设为常数. 使用单侧采集观测系统,共计模拟201炮和201道,炮间距和道间距分别为20 m和10 m. 表中的参数设置保证了各界面的速度扰动大小一致,能更好地验证角度域算法的准确度.

      表 1  常密度水平层状模型

      Table 1.  Horizontal layered model with constant density

      层厚度/m速度/(m·s–1)密度/(kg·m–3)
      30020002500
      30021002500
      30020002500
      30021002500
      30020002500
      30021002500

      图1a为基于水平层状模型数据计算的常规GRT反演剖面,界面位置准确,但浅层存在明显的过度叠加现象,这是因为不同位置的照明不均匀没有得到有效的处理. 图1b为角度域GRT反演剖面,可以看到,因为考虑了有效的照明范围,同相轴更加干净均匀. 在图1的两个剖面中,分别抽取三道(白色虚线位置)对比两种方法的反演振幅准确度. 图2显示了从常规GRT反演剖面中抽取的三道(黑色虚线)与真实扰动(红色实线)的对照图. 图3显示了从角度域GRT反演剖面中抽取的三道(黑色虚线)与真实扰动(红色实线)的对照图. 对比可见,角度域反演结果波形更加规整,振幅更加准确. 图4显示了两种方法的水平保幅情况,其中图4a是将图1a常规GRT反演剖面中水平界面的峰值振幅抽取出来,图4b是将图1b角度域GRT反演剖面中水平界面的峰值振幅抽取出来,可以看到,角度域方法获得了更宽的保幅范围.

      图  1  (a)常规GRT反演剖面;(b)角度域GRT反演剖面(修改自Li et al., 2018

      Figure 1.  Inversion profiles for the layered model by using (a) conventional GRT inversion method and (b) angle-domain GRT inversion method (modified from Li et al., 2018)

      图  2  常规GRT反演振幅曲线(黑色虚线)和真实扰动振幅(红色实线)对比图. 道集位置对应图1a中的白色虚线位置: (a)x=500 m;(b)x=800 m;以及(c)x=2000 m(修改自Li et al., 2018

      Figure 2.  Comparison of true perturbations (red solid lines) and retrieved perturbations (black dashed lines) by conventional GRT inversion method at the distances (a) x=500 m ; (b) x=800 m ; and (c) x=2000 m in Fig. 1a (modified from Li et al., 2018)

      图  3  与图2类似,角度域GRT反演振幅曲线与真实扰动振幅对比图(修改自Li et al., 2018

      Figure 3.  Similar to Fig. 2, but for angle-domain GRT inversion method from Fig. 1b (modified from Li et al., 2018)

      图  4  水平保幅情况.(a)常规GRT反演;(b)角度域GRT反演(修改自Li et al., 2018

      Figure 4.  The retrieved amplitude curves picked along the reconstructed reflectors from (a) conventional GRT inversion profile and (b) angle-domain GRT inversion profile (modified from Li et al., 2018)

      表2为水平多层双参数模型,其中各层密度与速度同步变化. 各界面的速度和密度扰动大小一致,能更好地验证算法的准确度. 模型大小网格参数以及地震采集参数与单参数测试基本一致.

      表 2  变密度水平层状模型

      Table 2.  Horizontal layered model with varying density

      层厚度/m速度/(m·s–1)密度/(kg·m–3)
      30020002000
      30021002100
      30020002000
      30021002100
      30020002000
      30021002100

      图5显示了由方程(18)计算的双参数水平层状模型的中间道集处角度域共成像点道集,为进一步观察生成的角度域共成像点道集的AVA特征和保幅特性,将角道集中的峰值曲线抽取出来并与理论曲线做对比. 图6a显示了横向距离x = 2 000 m处的AVA理论振幅曲线,横坐标是角度余弦值,纵坐标是振幅. 可以看到,5条曲线因为正负重叠成了2条曲线,这与界面参数扰动相对应. 另外在理论情况,假定各成像点拥有足够的照明覆盖,5条曲线呈线性平稳变化. 图6b则显示了角度域方法计算并经标定后的相同道集处的AVA振幅曲线. 可以看到在每个成像点有效的照明范围内,AVA特征曲线与理论值基本保持一致. 由此验证了角度域逆散射偏移算法的准确性. 与常密度情况有明显区别的是,这里的AVA偏移振幅曲线存在明显的角度依赖,在有效的照明范围内,偏移振幅随散射角的余弦值呈线性变化,其中截距和斜率分别与扰动参数f1f2有关.

      图  5  变密度水平层状模型角度域共成像点道集图(修改自Li et al., 2018

      Figure 5.  ADCIGs calculated at the distance x=2 000 m (modified from Li et al., 2018)

      图  6  (a)理论AVA特征曲线;(b)计算提取的AVA曲线(修改自Li et al., 2018

      Figure 6.  (a) Theoretical AVA curves and (b) Picked AVA curves from the calculated ADCIGs panel (Fig. 5) (modified from Li et al., 2018)

    • 前文回顾了散射场的Born近似和级数表示,经典的线性逆散射方法有意义地选择散射级数首阶项(一阶Born近似)来近似散射场,并在此基础上,利用逆GRT方法实现线性化逆散射反演或偏移. 一阶Born近似实际上是将总波场近似为背景介质的入射场,忽略了扰动引起的多次散射,该方法仅适用于小扰动情况(Weglein and Gray, 1983Bleistein,1987符力耘等,1998符力耘,2010),而且在复杂非均匀介质中,多次散射效应对散射场能量的贡献往往不能忽略. 为了克服单散射线性反演方法仅对弱散射介质有效的局限性,Mao等(2013)在线性逆散射GRT偏移反演方法的基础上提出一种基于二阶Born 近似的非线性反演新思路. 该方法尝试考虑散射序列高阶项来弥补单散射的不足,主要探讨散射场序列的前两项,即二阶Born近似:

      $\begin{split} &{U_{sc}} \left({{{r}},{{s}},{\rm{\omega }}} \right) \approx {\left({{U_s}} \right)_1} + {\left({{U_s}} \right)_2} = {{\rm{\omega }}^2}\int {{\rm{d}}{{x}}{G_0}\left({{{r}},{{x}},{\rm{\omega }}} \right)f\left({{x}} \right)}.\\& \left[ {{G_0}\left({{{x}},{{s}},{\rm{\omega }}} \right) + {{\rm{\omega }}^2}\int {{\rm{d}}{{y}}{G_0}\left({{{x}},{{y}},\omega } \right)f\left({{y}} \right)} {G_0}\left({{{y}},{{s}},{\rm{\omega }}} \right)} \right] \end{split}$

      (21)

      其中二次散射积分项为:

      $\begin{split} \left({{U_s}} \right){}_2 & = {{\rm{\omega }}^2}\int {{\rm{d}}{{x}}{G_0}\left({{{r}},{{x}},{\rm{\omega }}} \right)f\left({{ x}} \right){{\rm{\omega }}^2}}.\\&\qquad\int {{\rm{d}}{{y}}{G_0}\left({{{x}},{{y}},\omega } \right)f\left({{y}} \right)} {G_0}\left({{{y}},{{s}},{\rm{\omega }}} \right) \end{split} $

      (22)

      图7简要给出了以扰动点x为主散射点的Born序列示意图,s为炮点位置,r为检波点位置. 其中图7a为全波场散射路径图示,地震波分别与其他散射点经过不同阶次的散射后与主散射点x相互作用传播到r图7b为一阶Born近似散射路径图示,地震波从s出发,传播到x发生一次散射作用再传播到r,也称为单散射;图7c是二阶Born近似散射路径示意图,包括单散射(实线)和二次散射(虚线),其中二次散射是地震波与所有散射点发生一次散射作用,再经主散射点x散射一次,最后传播到r的过程.

      图  7  散射序列示意图.(a)全波场散射示意图.(b)一阶Born近似散射示意图.(c)二阶Born近似散射示意图(李武群等, 2017

      Figure 7.  Schematics of scattering series for (a) full order Born scattering, (b) first order scattering, and (c) second order scattering (Li et al., 2017)

      与基于单散射一阶Born近似的线性偏移反演方法相比,求解基于二阶Born近似的积分方程反问题的难点在于如何有效地处理二次散射项. Mao等(2013)提出局部二次散射的概念,假定二次散射的贡献主要集中在一次散射点附近的圆域中,如图8. 利用一阶泰勒近似,可以将圆形邻域${B_{{x}}}$内的散射点${{y}}$的散射位势表示为:

      图  8  局部二次散射示意图(李武群等,2017

      Figure 8.  Schematics of second order scattering within a local area (Li et al., 2017)

      $f\left({{y}} \right) \approx f\left({{x}} \right) + {\nabla _{{x}}}f\left({{x}} \right) \cdot \left({{{y}} - {{x}}} \right)$

      (23)

      此时,方程式(21)二阶Born近似${U_{sc}}$可以合理地近似表示成式(24)${\widetilde U_{sc}}$,根据$f\left({{x}} \right)$$c\left({{x}} \right)$的关系,式(24)可以进一步写为式(25),单独处理方括号中关于${{y}}$的积分部分可得式(26),该式是基于局部二次散射的二阶Born近似表达式,这个表达式类似于一阶Born近似(单散射)表达式. 推导出该式的目的是将二阶Born近似表达成类似声波GRT的结构,通过逆GRT偏移反演出关于散射位势$f$的二次项,进一步实现非线性二次反演. 式(26)中H为二次散射系数,只与背景模型速度以及局部二次散射半径因子$l$有关,可以看成是二次扰动的加权项,是局部二次散射假设的核心所在,具体推导过程可以参考Ouyang等(2014).

      ${\widetilde U_{sc}}\left({{{r}},{{s}},{\rm{\omega }}} \right) = {{\rm{\omega }}^2}\int {{\rm{d}}{{x}}{G_0}\left({{{r}},{{x}},{\rm{\omega }}} \right)f\left({{x}} \right){G_0}\left({{{x}},{{s}},{\rm{\omega }}} \right)} + {{\rm{\omega }}^2}\int {{\rm{d}}{{x}}{G_0}\left({{{r}},{{x}},{\rm{\omega }}} \right)f\left({{x}} \right) } {{\rm{\omega }}^2}\int_{{B_{{x}}}} {{\rm{d}}{{y}}{G_0}\left({{{x}},{{y}},\omega } \right)f\left({{y}} \right){G_0}\left({{{y}},{{s}},{\rm{\omega }}} \right)} $

      (24)

      ${\widetilde U_{sc}}\left({{{r}},{{s}},{\rm{\omega }}} \right) = {{\rm{\omega }}^2}\int {{\rm{d}}{{x}}{G_0}\left({{{r}},{{x}},{\rm{\omega }}} \right){G_0}\left({{{x}},{{s}},{\rm{\omega }}} \right)} \times \left[ {f\left({{x}} \right){\rm{ + }}{f^2}\left({{x}} \right){{\rm{\omega }}^2}\int_{{B_{{x}}}} {{\rm{d}}{{y}}{G_0}\left({{{x}},{{y}},\omega } \right)\left({1{\rm{ - }}\dfrac{{2{\nabla _{\bf{x}}}{c_0}\left({{x}} \right)}}{{{c_0}\left({{x}} \right)}} \cdot \left({{{y}} - {{x}}} \right)} \right)\dfrac{{{G_0}\left({{{y}},{{s}},{\rm{\omega }}} \right)}}{{{G_0}\left({{{x}},{{s}},{\rm{\omega }}} \right)}}} } \right]$

      (25)

      ${\widetilde U_{sc}}\left({{{r}},{{s}},{\rm{\omega }}} \right) \approx {{\rm{\omega }}^2}\int {{\rm{d}}{{x}}{G_0}\left({{{r}},{{x}},{\rm{\omega }}} \right){G_0}\left({{{x}},{{s}},{\rm{\omega }}} \right)\left[ {f\left({{x}} \right) + {f^2}\left({{x}} \right)H\left({{{x}},{{s}},{\rm{\omega }}} \right)} \right]} $

      (26)

      对方程式(26)应用逆GRT反演算子,可以得到在局部二阶Born近似意义下的关于扰动量$f$的二次多项式:

      $\left({{\Re ^ * }U} \right)\left({{z}} \right) \approx a\left({{z}} \right)\left[ {\left\langle {f\left({{z}} \right)} \right\rangle + H\left({{{z}},{{s}},l({{z}})} \right){{\left\langle {f\left({{z}} \right)} \right\rangle }^2}} \right]$

      (27)

      式中,z是成像点,${\Re ^ * }$是逆GRT偏移算子,$a\left({{z}} \right)$为散射夹角相关的积分. $U$为散射场数据. 特别地,在二维中,$U$为希尔伯特变换后的散射场数据. 通过该式,在线性反演的基础上,求解一个一元二次方程即可得到二阶Born近似意义下非线性解$\left\langle {f\left({{z}} \right)} \right\rangle $.

      散射场二阶Born近似考虑了二次散射的作用,理论上能模拟更高精度的地震波场. 图9为常密度单界面层状模型,水平方向和深度方向网格节点数分别为401和251,网格间距为10 m和5 m,扰动界面位置在深度1 000 m. 界面以上为常数背景速度${c_0}$=2 000 m/s,界面以下的速度依据速度扰动2%、10%、20%、30%、40%设置. 采用中间单炮激发,两侧接收,震源函数选择主频20 Hz的里克子波,采样时间4 s,采样间隔4 ms. 分别利用式(8)和式(21)计算这些不同扰动模型的一阶Born数据和二阶Born数据,并以解析数据为标准,分别计算这5个不同扰动模型的一阶、二阶Born数据振幅的相对误差. 图10为不同扰动模型的一阶Born数据与解析数据近偏移距道集振幅误差分布曲线图,图10b为不同扰动模型的二阶Born数据与解析数据近偏移距道集振幅误差分布曲线图. 由图10a可见,在小速度扰动(${c}$=2 040 m/s,δc=2%)时,一阶Born数据误差曲线接近0值并且随偏移距变化平缓,各道振幅误差都很小. 随着速度扰动由2%增大到40%,误差明显变大. 由图10b可见,不同扰动大小模型的二阶Born数据振幅误差都比较小,基本控制在5%以内,并且误差随偏移距变化不明显. 因此,一阶Born近似合成数据在小扰动情形是比较准确的,但随着速度扰动变大,误差相应地快速增大. 二阶Born近似数据准确度则明显提高,对于不同速度扰动的模型,合成数据都非常接近解析解.

      图  9  不同速度扰动的单界面水平层状模型

      Figure 9.  Single interface horizontal layered model with different velocity perturbations

      图  10  不同扰动模型Born近似误差变化曲线.(a)一阶Born近似误差.(b)二阶Born近似误差(李武群等,2017

      Figure 10.  Error distribution curves from Born data of different perturbation. (a) First order Born approximate data error. (b) Second order Born approximate data error (Li et al., 2017)

      图11为各扰动模型的线性、非线性反演剖面中界面对应的归一化振幅曲线. 由图11a可见,随着扰动量逐渐增大,GRT线性反演结果逐渐偏离真实值1(红色实线). 这表明基于一阶Born近似线性反演理论只适用于小扰动变化的介质模型. 图11b则是GRT非线性反演的振幅归一化曲线图. 由图可见,对于不同扰动程度的合成数据进行GRT非线性反演,都能获得比较准确的振幅.

      图  11  GRT线性反演和非线性反演界面振幅归一化曲线图.(a)GRT线性反演;(b)GRT非线性反演(李武群等,2017

      Figure 11.  Comparison of normalized inverted amplitudes from GRT linear and nonlinear inversion of different perturbation models. (a) GRT linear inversion. (b) GRT nonlinear inversion (Li et al., 2017)

      Ouyang等(2014)对横向变化复杂模型的有效性进行了测试. 图12a为Sigsbee 2A左侧部分. 图12b为GRT非线性反演剖面. 图13为水平坐标x=6 000 ft处(白色虚线)抽取的反演值与真实扰动值对比图. 其中图13a为GRT线性反演振幅曲线对比图,图13b为GRT非线性反演振幅曲线对比图. 实线表示反演值,虚线表示真实扰动值. 由图可见,在小扰动的区域(深度位于18 000 ft和29 000 ft之间),线性和非线性反演的振幅与真实值几乎一致;而在大扰动的位置(深度位于29 000 ft左右),非线性反演的结果更加准确.

      图  12  (a)Sigsbee 2A速度模型及(b)非线性反演剖面(修改自Ouyang et al., 2014

      Figure 12.  (a) the Sigsbee 2A model and (b) the nonlinear inversion profile (modified from Ouyang et al., 2014)

      图  13  反演值(实线)与真实扰动(虚线)之间的比较.(a)线性反演值与真实扰动之间的比较.(b)非线性反演值与真实扰动的比较(修改自Ouyang et al., 2014

      Figure 13.  Comparison between true perturbations (dotted line) and inverted perturbations (solid line). (a) Linear inversion. (b) Quadric nonlinear inversion (modified from Ouyang et al., 2014)

    • 地震勘探领域中涉及到逆散射级数方法的研究方向包括地球物理反演、地震深度成像、面波分析和多次波压制等. 其中多次波压制主要是利用逆散射级数的一些子序列来预测相应的自由表面多次波或层间多次波,然后在原始波场中减去预测多次波成分以实现多次波衰减压制. 本文从多次波压制方面介绍ISS法.

      根据Weglein等(2003),散射场可以由参考波场${G_0}$和扰动算子$V$以级数的形式给出,这里的扰动算子$V$可以表示成关于数据阶数的级数:

      $V{\rm{ = }}{V_1}{\rm{ + }}{V_2}{\rm{ + }}{V_3}{\rm{ + }} \cdot \cdot \cdot $

      (28)

      式中,${V_i}\left({i = 1,2,3,...} \right)$ 表示关于数据的第i阶成分. 若记$D$为散射场${U_s}$在观测表面上记录到的数据,将式(26)代入式(7)即可得到:

      $\begin{split} D =& {G_0}\left({\sum\limits_{i = 1}^\infty {{V_i}} } \right){G_0} + {G_0}\left({\sum\limits_{i = 1}^\infty {{V_i}} } \right){G_0}\left({\sum\limits_{i = 1}^\infty {{V_i}} } \right){G_0} + \\& \cdot \cdot \cdot + {G_0}{\left\{ {\left({\sum\limits_{i = 1}^\infty {{V_i}} } \right){G_0}} \right\}^n} + \cdot \cdot \cdot \end{split} $

      (29)

      ${V_{}}$级数展开,令等式两边各阶项相对应,有:

      $D = {G_0}{V_1}{G_0}$

      (30)

      $0 = {G_0}{V_2}{G_0}{\rm{ + }}{G_0}{V_1}{G_0}{V_1}{G_0}$

      (31)

      $\begin{split} 0 =& {G_0}{V_3}{G_0}{\rm{ + }}{G_0}{V_1}{G_0}{V_2}{G_0}{\rm{ + }}{G_0}{V_2}{G_0}{V_1}{G_0}{\rm{ + }}\\&{G_0}{V_1}{G_0}{V_1}{G_0}{V_1}{G_0} \cdots \end{split} $

      (32)

      根据以上方程组,利用散射场数据$D$求出${V_1}$,再由${V_1}$求出${V_2}$,以此类推,最终可以实现对逆散射级数${V_{}}$的求解. 以上推导过程对${G_{}}$${G_0}$并没有任何要求,也不涉及到小扰动的假设. 唯一的假定是,认为${V_1}$${V_{}}$的线性部分. 如果假设介质扰动很小,即${V_1}$接近${V_{}}$,那么根据式(30)即可求出${V_{}}$的近似解,此即前文中基于Born近似的线性化逆散射近似解.

      逆散射级数方法提供了一个基于参考介质直接求解多维非线性反演问题的框架,而且适应任意非均匀性. 以海洋采集为例,半空间水体对应的参考格林函数可表示为震源直接传播到接收点的波场$G_0^{\rm{d}}$和震源经过自由表面的反射后再传播到接收点的波场$G_0^{{\rm{fs}}}$之和:

      ${G_0}{\rm{ = }}G_0^{\rm{d}} + G_0^{{\rm{fs}}}$

      (33)

      假定$\overline D $为去掉检波点和炮点端鬼波的波场数据,与自由表面多次波有关的级数可以表示为:

      $D_1^{'} = \overline D = G_0^{\rm{d}}{V_1}G_0^{\rm{d}}$

      (34)

      $D_2^{'} = - G_0^{\rm{d}}{V_1}G_0^{{\rm{fs}}}{V_1}G_0^{\rm{d}}$

      (35)

      $\begin{split} D_3^{'} = &- G_0^{\rm{d}}{V_1}G_0^{{\rm{fs}}}{V_1}G_0^{{\rm{fs}}}VG_0^{\rm{d}} - G_0^{\rm{d}}{V_1}G_0^{{\rm{fs}}}{V_2}G_0^{\rm{d}} -\\& G_0^{\rm{d}}{V_2}G_0^{{\rm{fs}}}{V_1}G_0^{\rm{d}} \cdots \end{split} $

      (36)

      因此,去除全部自由表面多次波的数据为$D_{}^{'} = \displaystyle\sum\limits_{i = 1}^\infty {D_{{i}}^{'}},$其中还包括层间多次波和主反射波. 令$D_{}^{'}$中的${V_{}}$级数表示为:

      $V = \sum\limits_{i = 1}^\infty {V_{{i}}^{'}} $

      (37)

      进一步地,与层间多次波有关的级数可以表示为:

      $D_1^{''} = {D^{'}} = G_0^{\rm{d}}V_1^{'}G_0^{\rm{d}}$

      (38)

      $D_2^{''} = - G_0^{\rm{d}}V_1^{'}G_0^{{\rm{fs}}}V_1^{'}G_0^{\rm{d}}$

      (39)

      $\begin{split} D_3^{''} =& - G_0^{\rm{d}}V_1^{'}G_0^{{\rm{fs}}}V_1^{'}G_0^{{\rm{fs}}}V_1^{'}G_0^{\rm{d}} - G_0^{\rm{d}}V_1^{'}G_0^{{\rm{fs}}}V_2^{'}G_0^{\rm{d}} -\\& G_0^{\rm{d}}V_2^{'}G_0^{{\rm{fs}}}V_1^{'}G_0^{\rm{d}} \cdots \end{split} $

      (40)

      因此,去除全部自由表面多次波和层间多次波的数据为$D_{}^{''} = \displaystyle\sum\limits_{i = 1}^\infty {D_{{i}}^{''}},$即为余下的主反射波数据. 上述去多次波步骤比较灵活,先后顺序可以调换,也可以仅实现其中一种多次波的去除处理.

      ISS预测多次波方法是基于多维直接反演的思路,而且无需任何地下的先验信息. 理论上,该方法利用散射序列构建了地震波场与介质物性参数之间较为完整的非线性关系. 通过构建目标子序列能够准确预测目标多次波的走时和振幅信息. 当前工业标准的表面多次波去除方法(Surface-Related Multiple Elimination, SRME)首先构建自由表面多次波近似振幅和相位算子,然后利用能量最低自适应相减方法搭建预测多次波与实际多次波之间的联系. 该方法能够有效地去除受其他信号干扰小的表面多次波,如果目标表面多次波同相轴受到干扰,比如与反射波邻近,或者叠合的情形,能量最低自适应相减法会失效,而且会损伤邻近的有效反射波信号. 逆散射级数自由表面多次波衰减ISS-FSME(ISS Free-Surface Multiple Elimination)因为能够预测准确的走时和振幅,能够有效去除孤立或者受到干扰的表面多次波信号,最大程度保护有效反射能量. 针对层间多次波衰减问题,ISS方法通过应用能量最小化自适应去除孤立的目标层间多次波. 对于受到干扰的层间多次波,需结合更综合性的预测方法,确保邻近信号不受到影响.

      图14给出了墨西哥湾某区实际数据的ISS自由表面多次波去除例子,左图是原始数据的叠加剖面,右图为经ISS去除表面多次波后的叠加剖面,图15以偏移距域共成像点道集形式给出了墨西哥湾某区实际数据的ISS层间多次波去除例子. 中间线左侧显示了水平位置1 450 ft处的处理效果,中间线右侧显示了水平位置2 350 ft处的处理效果,分别从左到右为原始数据共成像点道集,预测的层间多次波成像点道集和去除层间多次波后的成像点道集.

      图  14  墨西哥湾某区ISS自由表面多次波去除应用示例(修改自Weglein et al., 2003

      Figure 14.  The left panel is a stack of a field data set from the Gulf of Mexico. The right panel is the result of inverse-scattering free-surface multiple removal (modified from Weglein et al., 2003)

      图  15  墨西哥湾某区ISS层间多次波去除应用示例(修改自Weglein et al., 2003

      Figure 15.  An example of inverse-scattering internal multiple attenuation from the Gulf of Mexico (modified from Weglein et al., 2003)

      逆散射级数法的优势在于不依赖地下介质信息也能预测多次波运动学和动力学信息,在现代多次波处理中受到高度关注. 散射序列是一种无穷序列,目标子序列的构建和选取需要精巧的设计,如何选取其中的子序列应用于实践中相应的问题,需要更多的思考. 高维复杂介质数值计算的稳定性以及计算成本也是影响其进一步推广应用的重要原因.

    • 当前主流的地震偏移成像技术有Kirchhoff偏移,高斯束偏移和逆时偏移. Kirchhoff偏移成像方法具有灵活、高效、低门槛等特点,是工业界普遍采用的偏移成像方法. 高斯束偏移成像(GBM)方法以能解决多路径、多走时问题和陡倾界面成像的能力著称,可以有效对复杂构造焦散区及盲区成像. 逆时偏移成像技术则是当前公认的精度最高的成像方法,具有利用回转波进行成像的优点. 基于逆散射反演理论的成像方法,比如GRT偏移成像方法,属于基于射线理论的波动方程积分解法,拥有射线类成像的共同优点. 从某种意义上,GRT偏移可以视为广义Kirchhoff偏移的一个加权版本(Dillon, 1990),具有优于Kirchhoff偏移的保幅特性和优于高斯束偏移和逆时偏移的计算效率. 当前,关于逆散射偏移成像(ISM)的研究和应用相对较少,本节通过一些合成数据和实际资料的测试和应用算例介绍ISM成像方法的成像效果和特点.

      (1)合成数据测试

      图16为SEAM工业标准模型东西向的一个2D速度剖面,该模型中间为一高速盐丘,盐翼较陡倾,两侧为较平缓的沉积层. 模型水平和深度方向坐标范围分别为0~35 km和0~15 km. 声波波速分布如图,依照Gardner关系式$\rho {\rm{ = }}0.31 \cdot {V_p}_{}^{0.25}$(波速单位:m/s,密度单位:g/cm3)设定模型密度. 选择交错网格有限差分方法计算合成数据(Thorbecke and Draganov, 2011),为了简便起见,原模型网格间隔被离散成10 m. 炮点检波点布置在海面,在0~35 km范围内布置共计876炮以及3 501个检波点,炮间距和道间距分别为40 m和10 m. 正演计算采用最大频率25 Hz的里克子波,记录时长8 s.

      图  16  二维SEAM模型

      Figure 16.  The 2D SEAM model

      图17为二维SEAM模型的偏移成像结果,其中图17ab分别为逆散射偏移成像(ISM)和高斯束偏移成像(GBM). 两个偏移剖面经过标定后使用相同的振幅范围. 由图可见,逆散射偏移(ISM)剖面振幅整体更加均衡,盐腹下方地层连续性更好.

      图  17  2D SEAM 模型偏移成像结果.(a)逆散射成像;(b)高斯束成像

      Figure 17.  The migration results of the 2D SEAM model by (a) ISM, and (b) GBM

      为了验证逆散射偏移成像方法在3D地震数据中的适用性,本文以3D SEG/EAGE盐丘模型为例进行测试,模型如图18. 该模型是常密度声波介质模型,整个区域内薄砂层、断层遍布,中间区域包含一高速盐丘,这些地质构造叠加在一个光滑背景速度波形上. 模型东西横向和南北纵向坐标范围均为0~13.5 km,深度范围0~4.2 km,模型沿横向、纵向和深度方向的网格点数分别为676、676和210,网格间隔均为20 m.

      图  18  3D SEG/EAGE盐丘模型

      Figure 18.  The 3D SEG/EAGE salt model

      为方便说明,本文抽取三维偏移体中的几个切片来分析偏移结果,分别选择横坐标方向和纵坐标方向共计4个竖向剖面以及1个水平深度切片,并将偏移结果与高斯束偏移成像比较,结果如图所示. 图19显示了纵坐标方向两个剖面的成像结果,其中图19ab分别为y=7.8 km剖面和y=6.6 km剖面的原始速度模型. 图19cd为对应的逆散射偏移(ISM)剖面,图19ef为对应的高斯束偏移(GBM)剖面. 由图可见,对于主要的盐丘体结构,薄砂层以及断层等构造,两种偏移方法均获得了较好的成像. 逆散射偏移(ISM)对于结构的刻画更加精细,分辨率更高,高斯束偏移(GBM)则拥有更好的信噪比. 图20显示了横坐标方向两个剖面的成像结果,其中图20ab分别为x=5.28 km剖面和x=2.04 km剖面的原始速度模型. 图20cd为对应的逆散射偏移(ISM)剖面,图20ef为对应的高斯束偏移(GBM)剖面. 与横坐标方向的剖面结果类似,两种偏移方法成像效果基本相当. 逆散射方法对于盐丘体(边界、内部),较浅层结构等区域成像保幅一致性较好,分辨率稍高. 高斯束方法有考虑多路径的优势,整个成像剖面信噪比更好. 图21a为深度方向z=0.8 km处的原始速度切片. 图21b为逆散射偏移(ISM)切片,图21c为对应的高斯束偏移(GBM)切片. 从成像效果来看,两种偏移方法差异并不明显,逆散射偏移(ISM)切片分辨率稍高. 由于数据的限制原因,周边区域无法成像.

      图  19  3D SEG/EAGE模型沿x方向剖面偏移成像结果.(a)和(b)分别原始速度模型.(c)和(d)为逆散射偏移成像结果.(e)和(f)为高斯束偏移成像结果

      Figure 19.  Migration results of the profiles along x direction. (a) and (b) are velocity models. (c) and (d) are ISM results. (e) and (f) are GBM results

      图  20  3D SEG/EAGE模型沿y方向剖面偏移成像结果.(a)和(b)分别原始速度模型.(c)和(d)为逆散射偏移成像结果.(e)和(f)为高斯束偏移成像结果

      Figure 20.  Migration results of the profiles along y direction. (a) and (b) are velocity models. (c) and (d) are ISM results. (e) and (f) are GBM results

      图  21  3D SEG/EAGE模型深度切片偏移成像结果. (a)为深度位置z=0.8 km切片上的速度模型. (b)为逆散射偏移成像结果. (c)为高斯束偏移成像结果

      Figure 21.  Migration results of the depth slice. (a) Velocity at slice z=0.8 km. (b) ISM result. (c) GBM result

      (2)实际资料测试

      本文对某工区的实际采集数据进行了逆散射成像处理. 图22是该工区长排列采集观测系统俯视图,图中的蓝色小三角形表示检波点,红色小圆点表示炮点. 该区采集布置了一条炮线和两条检波线,检波线间距距离不超过200 m,其中有少许炮点与检波点非常零散地分布在检波线附近. 所有炮点与检波点在纵向距离为300 m的范围内,测线全长32 km. 共计激发828炮,每炮激发所有检波点接收,共3 202个检波点,总计2 651 256道,每道地震记录有2 500个采样点,采样间隔为4 ms. 图23为该工区野外采集的某单炮地震记录,可以看到,记录中有效反射同相轴不明显,噪声较大. 图24为逆散射偏移成像剖面和局部放大图,图25为Kirchhoff积分偏移成像剖面和局部放大图. 对比虚线方框可以看到,逆散射成像对断层细节刻画得更突出. 另外,对于深部构造的呈现也更加明显.

      图  22  长排列采集观测系统

      Figure 22.  The Schematic of the acquisition geometry

      图  23  野外单炮地震记录

      Figure 23.  The single shot record of the field data

      图  24  逆散射偏移成像剖面(左)和局部放大区域(右)

      Figure 24.  The imaging result by Inverse scattering migration (left) and the zoomed one (right)

      图  25  Kirchhoff偏移成像剖面(左)和局部放大区域(右)

      Figure 25.  The imaging result by Kirchhoff migration (left) and the zoomed one (right)

      (3)VSP数据逆散射成像

      垂直地震剖面(Veritcal Seismic Profile, VSP)观测系统是一种在地表激发震源,检波点置于井孔中进行接收的地震采集方法. 相比地面地震资料,井中采集避开了浅层低速带的高频衰减和地表的噪声干扰,获取的资料具有更高的信噪比和分辨率,接收到的波场更加丰富,而且能够直接获取更具保真度的地层参数. 这些固有的特性使得VSP技术在获得地下井旁构造的高精度高分辨率成像方面更具潜力. Miller等(1987)最早测试了逆散射GRT偏移方法对VSP数据的真振幅重建,其给出的算例基于均匀背景速度,偏移中也只用到了上行反射波信号. 以此为基础的常规逆散射VSP成像方法也主要利用上行反射波信号进行成像. 当前主流的VSP成像方法,尤其生产中使用的方法都是基于上行反射波信号成像的方法. 这些方法广泛存在波场分离不彻底,滤波参数不好控制等问题. 另外,仅利用上行反射波信息,往往只能对井旁局部范围照明,大量多次波有效信号没有得到合理利用. 新近出现了一些利用下行多次波进行成像的方法,实现了更宽范围的成像. 其中射线类偏移方法仍需要在波场分离后再对多次波信号成像,而且单一强调多次波的效果,忽略上行反射波和下行多次波信号的优缺点和互补的重要性.

      基于VSP波场信号照明特点和上下行波场干涉特征的认识,Li等(2019)提出了基于逆散射偏移的VSP全波成像方法. 该技术可以直接对全波场VSP数据进行处理,无需进行波场分离,在提高效率的同时避免因波场分离带来的问题. 该方法可以联合上、下行波信号进行全波成像,发挥VSP资料的最大优势,改善常规VSP成像技术的井旁照明效果,实现井旁及远井区域的宽范围成像.

      图26a为SEAM P波速度模型的一个二维切片. 该二维速度切片取自SEAM模型北向23 900 m处. 左侧区域为较平缓的沉积层,右侧红色区域为一高速盐丘,盐翼严重陡倾. 模型水平长度7 000 m,深度5 000 m,模型网格数701×501,网格间距10 m×5 m. 炮点布置在地表,在0~7 000 m范围内按均匀间隔设置176炮,炮间距40 m. 垂直井位置在距离为3 000 m处(白色竖线). 井深范围0~4 000 m,检波器设置在井上,在0~4 000 m范围内按均匀间隔设置401道检波器,道间距10 m. 震源选择里克子波模拟,最大频率50 Hz. 记录长度8 s,时间采样间隔4.2 ms. 合成数据使用有限差分交错网格程序计算,模型顶部采用自由表面边界条件,四周和底部采用吸收边界.

      图  26  SEAM成像结果. (a),(b),(c)和(d)依次为速度场,上行波成像,下行波成像以及全波成像结果(修改自Li et al., 2019

      Figure 26.  Imaging results of the SEAM 2D model (a). (b), (c) and (d) are the imaging results using primary, multiple and joint imaging condition, respectively (modified from Li et al., 2019)

      该方法直接将SEAM 2D VSP全波场数据(含直达波)作为输入,分别给出上行反射波偏移,下行多次波偏移以及联合全波偏移的成像结果. 图26b为常规上行反射波偏移结果,可以看到井旁沉积层结构得到了清晰准确成像,包括右侧盐翼也得到了较好成像. 图26c为下行多次波偏移成像结果,对比图26bc虚线方框位置可以看到井旁成像范围得到明显拓宽. 另外对比图26bc白色箭头标记处,可以看到上行反射成像对近井盐翼结构的重建效果更好. 图26d是联合上下行波场全波成像结果,由于同时考虑了两种信号的优点,获得了最佳的成像效果.

    • 本文回顾了地震波散射和逆散射方法理论的发展和基本原理. 从经典线性逆散射GRT偏移反演技术出发,分别介绍了角度域逆散射反演方法,基于二阶Born近似和散射级数的非线性方法,以及基于逆散射理论的VSP成像方法等. 通过具体例子分别介绍了逆散射技术在合成数据和实际资料中的典型应用和最新进展.

      地震逆散射理论是建立在扰动理论和高频射线理论基础上的包含地震波场正演模拟和反演问题求解的理论体系. 具有严密的数学理论基础和生产应用上的灵活高效优势,而且可以描述勘探领域大部分理论和应用问题,并且找到相应的求解办法. 从成像效果来看,与当前主流的Kirchhoff、高斯束成像的射线类成像对比并不逊色,甚至具有更好的局部分辨率. 另外在线性还是非线性问题描述上,都有一套完备的反演理论体系,具有很大的挖掘空间应用前景. 结合逆散射方法的特点以及当前地震勘探发展的形势,笔者认为可以在以下方面进一步发挥逆散射方法的优势特点:

      (1)利用角度域逆散射偏移方法进行偏移速度分析,优化速度建模效率.

      (2)发展三维宽方位角度域保幅成像或反演方法,为五维地震资料处理和解释提供支持.

      (3)解决多波多分量弹性波地震处理中的问题,实现弹性波逆散射偏移成像方法的推广应用.

      (4)优化逆散射级数预测多次波方法,进一步提高多次波预测效率.

      (5)挖掘散射理论对于小尺度构造刻画的能力,发展高精度绕射波成像.

      笔者认为随着方法研究的推进和计算能力的提升,地震逆散射偏移与反演技术将逐步完善,实现从线性到非线性,从声波、弹性波到各向异性,从常规三维到五维,从常规地面采集到非常规采集(VSP/OBC/OBN)的研究延伸. 在此基础上,逆散射方法有望得到更广泛的应用.

参考文献 (49)

目录

    /

    返回文章
    返回