引用本文
邓文志, 李振春, 王延光, 孙小东. 基于稳定逆时传播算子的黏声介质最小二乘逆时偏移[J]. 物探与化探, 2015,39(4): 791-796.
DENG Wen-Zhi, LI Zhen-Chun, WANG Yan-Guang, SUN Xiao-Dong. The least-squares reverse time migration for visco-acoustic medium based on a stable reverse-time propagator[J]. Geophysical and Geochemical Exploration, 2015,39(4): 791-796.  
Doi:10.11720/wtyht.2015.4.22
Permissions
Copyright©2015, 物探与化探编辑部
? 《物探与化探》编辑部
基于稳定逆时传播算子的黏声介质最小二乘逆时偏移
邓文志1, 李振春1, 王延光2, 孙小东1
1.中国石油大学(华东) 地震波传播与成像实验室,山东 青岛 266580
2.中国石化胜利油田分公司 物探研究院,山东 东营 257068

作者简介: 邓文志(1989-),男,湖南人,现为中国石油大学(华东)硕士研究生,主要从事黏介质地震波传播与成像方面的研究工作。

收稿日期: 2014-09-10
基金: 国家自然科学基金项目(41374122); 国家科技重大专项项目(2011ZX05006-002)
摘要

基于GSLS模型黏声介质二阶拟微分方程,采用伪谱法进行数值模拟。针对黏声介质逆时传播过程中产生的高频不稳定问题,提出加入规则化算子对其进行消除的方法,构建了稳定的逆时传播算子。在最小二乘反演的基础上,将黏声介质逆时偏移与最小二乘思路相结合,发展了带有振幅补偿的黏声介质最小二乘逆时偏移(LSRTM)。Marmousi模型结果表明:相对于常规最小二乘逆时偏移,黏声介质最小二乘逆时偏移校正了地层的黏滞性,得到了更加精确可靠的保幅成像剖面。

关键词: 黏声介质; 拟微分方程; 逆时传播算子; 最小二乘逆时偏移
中图分类号:P631.4 文献标志码:A 文章编号:1000-8918(2015)04-0791-06
The least-squares reverse time migration for visco-acoustic medium based on a stable reverse-time propagator
DENG Wen-Zhi1, LI Zhen-Chun1, WANG Yan-Guang2, SUN Xiao-Dong1
1.SWPI,China University of petroleum (East China),Qingdao 266580,China
2.GRI of SINOPEC Shengli Oilfield Branch,Dongying 257068,China
Abstract

Based on a single SLS isotropic medium second-order pseudo-differential equation,the authors used the pseudo-spectral method to calculate numerical simulation.The authors propose introduction of regularization operator method to eliminate high-frequency instability caused by the backward propagation inverse process so as to build a stable reverse-time propagator and to achieve the amplitude compensating visco-acoustic medium least-squares reverse time migration.Marmousi model shows that,compared with least-squares reverse-time migration visco-acoustic least-squares reverse-time migration can correct the effect of viscosity and produce more accurate and better reliable amplitude preservation imaging section.

Keyword: visco-acoustic medium; pseudo-differential equation; reverse-time propagator; least-squares reverse-time migration

随着地震勘探日趋向着深层复杂地区发展, 面向复杂地区成像的地震处理成像方法研究越来越受到重视。由于地球广泛存在波动黏滞性[1, 2], 黏滞性呈现为地震波振幅随传播距离的增加而呈指数规律衰减[3]。为了获得高精度保幅成像, 需要校正黏滞性对成像的影响, 因此, 用黏滞性介质表示实际地层, 通过数值模拟方法来研究地震波在其中的吸收衰减规律, 进而指导地震资料处理就变得尤为重要。

黏介质的偏移成像是随着黏介质模型基本理论的发展而发展的, 最初人们应用反Q滤波技术对黏滞性进行补偿[3, 4, 5]; 随后, Dai等[6]实现了二维黏声介质叠后深度偏移; Traynin等[7]和Xie等[8]基于射线理论, 发展了叠前克希霍夫反Q偏移。由于在频率域很方便地进行衰减补偿, 国内外很多学者对黏介质单程波偏移进行了深入的研究[9, 10, 11, 12]。逆时偏移技术[13]作为一种精确的偏移方法成为现在主流偏移手法之一, 近年来许多学者将逆时偏移扩展到黏声介质中[14, 15, 16, 17]。Deng等[14]基于达朗贝尔黏声介质进行了保幅叠前深度偏移研究; Fletcher等[15]运用声波方程进行波场延拓, 在成像之前进行振幅和相位校正, 实现了黏声介质逆时偏移; Zhang等[16]基于线性黏声介质的频散关系推导出拟微分方程, 为了解决逆时偏移逆时传播过程中产生的高频不稳定问题, 考虑引入规则化算子, 不过并没有详细说明; Zhu等[17]修改拟微分方程并利用高频滤波器解决不稳定性问题, 也得到了较好的效果。

传统的成像理论并不是地震波正演算子的逆, 仅仅是它的共轭转置[18], 所以得到的只是近似解, 不能满足人们对复杂目标区域勘探开发的要求。最小二乘反演[19, 20, 21]能够帮助地震数据成像减少噪声、提高信噪比、平衡振幅, 得到比传统成像方法质量更高的成像结果。近年来, 很多学者从反演角度出发, 将逆时偏移与最小二乘思路相结合, 发展了在最小二乘反演框架下的逆时偏移(LSRTM)成像理论方法[22, 23, 24, 25], 相对于逆时偏移成像, 最小二乘逆时偏移成像可以克服RTM对深部储层成像分辨率不高且振幅不均衡等问题。为了克服最小二乘逆时偏移耗时太大的问题, 一些学者将相位编码技术[25, 26, 27, 28, 29]引入最小二乘逆时偏移中, 通过编码超道集的形成, 在保证成像精度的同时, 大大减少了计算耗时。

在以上基础上, 首先基于GSLS模型黏声介质二阶拟微分方程, 笔者采用伪谱法进行数值模拟, 推导了其伪谱法差分格式。然后针对黏声介质逆时传播过程中产生的高频不稳定问题, 提出加入规则化算子对其进行消除的方法, 构建了稳定的逆时传播算子。最后在最小二乘反演的基础上, 将黏声介质逆时偏移与最小二乘思路相结合, 发展了带有振幅补偿的黏声介质最小二乘逆时偏移成像。

1 方法原理
1.1 黏声介质拟微分方程

基于GSLS模型的二维黏声介质拟微分方程[30]

∂2∂t2+τv-Δ2∂∂t-v2ΔP=0。(1)

其中: 为拉普拉斯算子; v为传播速度; τ σ 和τ ε 分别代表应力和应变松弛时间, τ =τ ε /τ σ -1。应力和应变松弛时间可以表示为[31]

τσ=Q2+1-1wQ, τs=Q2+1+1wQ。(2)

可以看到, 拟微分方程相对于完全理想弹性介质声波方程而言, 多了一项衰减项, 当衰减足够小时, 可以退化到声波方程。

笔者运用伪谱法进行数值模拟, 伪谱法是时间导数采用有限差分求解, 空间导数利用傅式微分性质求解[32]。目前, 扩展的伪谱法将傅氏微分性质扩展, 可以用于求解空间分数阶偏导数[33]。对于一个波函数u(x), 通过傅式微分性质, 则其一阶、二阶偏导数可表示为

∂xu(x)=F-1[ikxF(u)], (3)∂x2u(x)=F-1(ikx)2F(u)]。(4)

对傅式微分性质进行扩展, 则空间分数偏导数表达式为

∂xβu(x)=F-1(ikx)βF(u)], (5)

这里F、F-1分别为傅里叶正变换、逆变换, β 为实数。

利用伪谱法对式(1)进行求解

p(t+2Δt)=2p(t+Δt)-p(t)+dt2{v2F-1[(ikx)2+(ikz)2]F[p(t+Δt)]}-τv2F-1{[-(ikx)2-(ikz)2]12F[p(t+Δt)-p(t)dt]}}。(6)

式(6)就构成了求解黏声介质伪谱法差分格式。

1.2 稳定的逆时传播算子构建

黏声VTI介质地震波传播是一个振幅能量衰减的过程, 为了进行偏移成像, 需要在接受波场逆时外推时, 补偿从反射界面到检波器传播过程中产生的振幅损失。因此, 需要修改正演方程式(1), 使方程衰减项符号变号[17], 此时, 逆时传播方程变为

∂2∂t2-τv-Δ2∂∂t-v2ΔP=0。(7)

比较黏介质单程波补偿算子[9]和基于旅行时的衰减补偿算子[14]可知, 基于双程波的算子能更方便地对地震波传播中引起的能量损失进行补偿, 还能体现逆时偏移对于高陡构造成像的优势。

然而, 逆时传播时会产生额外的高频成分, 呈现指数增加, 导致传播不稳定, 且Q值越小, 越不稳定。为了避免高频不稳定, 笔者借鉴Liu[34]的思想, 引入规则化算子( σv∂∂tΔ P)解决高波数引起的有限差分数值频散, 对高频不稳定问题进行消除。此时稳定的逆时传播算子为

∂2∂t2-τv-Δ2∂∂t-v2Δ+σv∂∂tΔP=0　。(8)

1.3 黏声介质最小二乘逆时偏移

对于黏声介质最小二乘逆时偏移, 基于Born近似下的线性化正演算子(即反偏移算子)可以写成

∂2∂t2+τs0-Δ2∂∂tp0(x, t)=f(t), ∂2∂t2+τs0-Δ2∂∂t-s02Δp(x, t)=m(x)∂2p0(x, t)∂t2。9

式中:s0为背景慢度模型, m(x)为扰动慢度模型, f(t)为初始点震源, p0(x, t)为背景波场, p(x, t)为扰动波场。扰动波场p(x, t)是通过两个过程计算得到:一是初始点震源f(t)和背景慢度s0 产生p0(x, t); 第二个也用到了背景慢度模型s0, 但是震源项是m(x)p0(x, t)/2t。

为了表示方便, 用矩阵向量的形式表示算子, 非线性正演算子定义为A , 线性化正演算子定义为L, 逆时偏移算子是LT。线性化正演算子L是A的Frechet 导数, LT是L的共轭。此时, 式(1)表示的正演过程为d=Am, 式(9)表示的反偏移过程为d=Lm。最小二乘反演的目的是找到一个扰动慢度m, 使得输入数据d=Lm满足误差函数最小

f(m)=12‖Lm-d‖2, (10)

其迭代解[23]为

m(k+1)=m(k)-αg(k), g(k)=LT[L(m(k))-d], α=(g(k))Tg(k)(Lg(k))T(Lg(k))。　　　　　　　　　(11)

这里, α 为迭代步长, g为梯度。迭代解可以用两种预处理算子进行约束, 第一种是空间域的高通滤波器[35]。在最小二乘迭代中, 需要定义一个合适的波数域的滤过带(通带), 以便区分假构造和模型更新量。然后在前几次迭代中应用高通滤波器来去除低波数背景散射假象。注意的是, 迭代共轭梯度算法中, 预处理算子必须是对称正定(SPD)的。然而, 不能自动化地确保用到的高通滤波器是SPD 的, 因此这种预处理算子在5次迭代之后一般便不再采用。第二种预处理算子是照明补偿, 即对角Hessian逆的一种近似[36]。这种照明补偿是一个元素全为正的对角矩阵, 因此它是对称正定的, 而且在每次迭代中都能放心使用, 因此, 文中采用照明补偿来进行预处理。我们可以将上述过程描述为图1所示流程。

图1 黏介质最小二乘逆时偏移流程

2 模型试算

首先对一个均匀黏声介质模型进行模型试算。模型纵、横向采样点数为400× 400, 纵、横向网格间隔为dx=dz=5 m, 速度参数为vp=2 000 m/s, 品质因子分别为无穷(完全理想弹性介质)、50、20, 震源为主频20 Hz的子波, 时间采样间隔为0.5 ms。图2为不同Q值下的波场快照及频谱图, 从图2可见:地层的黏滞性会损耗地震波能量, 使得主频降低, 频带宽度变窄, 说明本文的拟微分方程能准确地描述黏介质中地震波的运动学特征。

图2 不同Q值下的波场快照及频谱

接下来对逆时传播稳定性进行分析。图3是Q=50时(300 m, 300 m)处正向和逆时传播波形图及频谱图。由图3a和图3b可知, 当Q较小时, 逆时传播时会导致传播不稳定, 这种不稳定主要产生在高频成分。由图3c、图3d可知, 加入规则化算子后, 逆时传播高频不稳定现象得到消除。

图3 Q=50时(300 m, 300 m)处规则化前后正向和逆时传播的波形及频谱图

下面对Marmousi模型进行试算, Marmousi模型的P波速度场和Q模型如图4所示, P波速度范围为1 500~5 500 m/s, Q模型根据速度场建构, Q值范围为50~500, 在速度为2 000~4 000 m/s范围内存在强的衰减, 模型纵、横向采样点数为737× 380, 纵、横向网格间隔为dx=dz=10 m, 震源为主频20 Hz的子波, 总共184炮, 炮间距为40 m, 跑点范围是从第一个网格点到736个网格点。每炮检波器个数为737, 检波器间隔10 m, 时间采样为6.0 s, 时间采样间隔为1.0 ms, 炮记录采用伪谱法正演模拟得到。

通过对逆时传播加入规则化算子, 构建稳定的逆时传播算子, 进行黏声介质最小二乘逆时偏移成像。为了看衰减的影响, 对相同的黏声介质炮记录进行成像, 黏声介质LSRTM采用黏声波方程, 常规LSRTM使用声波方程, 都采用震源照明补偿作为预处理算子约束。图5为Marmousi模型偏移结果, 从图中可以看出, 常规LSRTM和黏声LSRTM都使得表层层位清晰可见, 地下轮廓明显, 构造形态准确, 没有低频噪声的影响, 在没有衰减的浅层, 有着相同的结果。然而在存在强的衰减的深层, 黏声最小二乘偏移由于补偿了反射界面到检波器传播过程中产生的能量损失, 相对于常规最小二乘偏移, 深部能量强, 振幅分布更加均衡, 得到更为精确可靠的保幅成像剖面。

图4 Marmousi模型参数

图5 Marmousi模型偏移结果

3 结论

基于GSLS模型黏声介质二阶拟微分方程, 采用伪谱法进行数值模拟, 针对黏声介质逆时传播过程中产生的高频不稳定问题, 提出加入规则化算子对其进行消除的方法, 构建了稳定的逆时传播算子。将基于双程波动方程数值解法的逆时偏移引入最小二乘偏移思想框架下, 充分结合了逆时偏移具有无倾角限制, 可以对转换波、多次波、棱柱波等特殊波进行成像, 而且成像方法不受介质速度变化的影响, 具有能很好处理纵、横向变速问题的优点与最小二乘偏移能降低成像假象、均衡反射幅值、提高分辨率的优点, 发展了带有振幅补偿的黏声介质最小二乘逆时偏移成像。结果表明, 黏声介质最小二乘逆时偏移补偿了地震波的吸收, 成像剖面整体效果清晰, 能量均衡, 深部能量得到加强, 使成像目标更接近于地质实际, 提高了深层构造的精细保真处理能力。

The authors have declared that no competing interests exist.

参考文献
[1] Carcione J M. Wave propagation in anisotropic linear viscoelastic media: Theory and simulated wavefields[J]. Geophysical Journal International, 1990, 101(3): 739-750. [本文引用:1]
[2] 孙成禹. 地震波传播理论与方法[M]. 东营: 中国石油大学出版社, 2007. [本文引用:1]
[3] Bickel S H, Natarajan R R. Plane-wave Q deconvolution[J]. Geophysics, 1985, 50(6): 1426-1439. [本文引用:2]
[4] Hargreaves N D, Calvert A J. Inverse Q filtering by Fourier transform[J]. Geophysics, 1991, 56(4): 519-527. [本文引用:1]
[5] 王珊, 于承业, 王云专, 等. 稳定有效的反Q滤波方法[J]. 物探与化探, 2009, 33(6): 696-699. [本文引用:1]
[6] Dai N, West G F. Inverse Q migration[C]//Expand ed Abstracts of the 64th SEG Annual International Meeting, 1994: 1418-1421. [本文引用:1]
[7] Traynin P J, Reilly J M. Amplitude and band width recovery beneath gas zones using Kirchhoff prestack depth Q-migration[C]//Expand ed Abstracts of the 78th SEG Annual International Meeting, 2008: 2412-2416. [本文引用:1]
[8] Xie Y, Xin K, Sun J, et al. 3D prestack depth migration with compensation forf requency dependent absorption and dispersion[C]//Expand ed Abstracts of the 79th SEG Annual International Meeting, 2009: 2919-2922. [本文引用:1]
[9] Zhang J, Wu J, Li X. Compensation for absorption and dispersion in prestack migration: An effective Q approach[J]. Geophysics, 2013, 78(1): S1-S14. [本文引用:2]
[10] 杨午阳. F-X 域黏弹性波动方程保幅偏移[J]. 石油物探, 2003, 42(3): 285-288. [本文引用:1]
[11] 孙天真, 谷玉田, 张惠欣, 等. 基于黏声介质的反Q滤波叠前深度偏移方法研究[J]. 石油物探, 2013, 52(3): 275-279. [本文引用:1]
[12] 郭恺, 娄婷婷. 双复杂介质条件下的反Q滤波偏移延拓算子研究[J]. 物探与化探, 2014, 38(3): 571-576. [本文引用:1]
[13] 屈念念, 李家斌. 基于P+S 波震源的弹性波叠前逆时偏移方法[J]. 物探与化探, 2013, 37(5): 859-865. [本文引用:1]
[14] Deng F, McMechan G A. True-amplitude prestack depth migration[J]. Geophysics, 2007, 72(3): S155-S166. [本文引用:3]
[15] Fletcher R P, Nichols D, Cavalca M. Wavepath-consistent effective Q estimation for Q compensated reserve-time migration[C]//Expand ed Abstracts of the 74th EAGE Annual International Conference and Exhibition, 2012. [本文引用:2]
[16] Zhang Y, Zhang H. Compensating for visco-acoustic effects in reverse-time migration[C]//Expand ed Abstracts of the 80th SEG Annual International Meeting, 2010: 3160-3164. [本文引用:2]
[17] Zhu T Z, Harris J M, Biondi B. Q-compensated reverse-time migration[J]. Geophysics, 2014, 73(9): S77-S87. [本文引用:3]
[18] Claerbout J F. Earth soundings analysis: Processing versus inversion[M]. Blackwell Science, 1992. [本文引用:1]
[19] Bamberger G, Chavent C, Hemon, et al. Inversion of normal incidence seismograms[J]. Geophysics, 1982, 47(5): 757-770. [本文引用:1]
[20] 杨其强, 张叔伦. 最小二乘傅立叶有限差分偏移[J]. 地球物理学进展, 2008, 23(2): 433-437. [本文引用:1]
[21] Chavent G, Plessix R E. An optimal true-amplitude least-squares pre-stack depth-migration operator[J]. Geophysic, 1999, 64(2): 508-515. [本文引用:1]
[22] Yao G, Jakubowicz H. Non-linear least-squares reverse-time migration[C]//Expand ed Abstracts of the 82nd SEG Annual International Meeting, 2012. [本文引用:1]
[23] Dai W, Schuster G T. Plane-wave least-squares reverse time migration[J]. Geophysics, 2013, 78(4): S165-S177. [本文引用:2]
[24] Zhang D L, Schuster G T. Least-squares reverse time migration of multiples[J]. Geophysics, 2014, 79(1): S11-S21. [本文引用:1]
[25] 黄建平, 曹晓莉, 李振春, 等. 最小二乘逆时偏移在近地表高精度成像中的应用[J]. 石油地球物理勘探, 2014, 49(1): 107-112. [本文引用:2]
[26] Dai W, Huang Y S, Schuster G T. Least-squares reverse time migration of marine data with frequency-selection encoding[J]. Geophysics, 2013, 78(4): S233-S242. [本文引用:1]
[27] Dai W, Wang X, Schuster G T. Least-squares migration of multisource data with a deblurring filter[J]. Geophysics, 2011, 76(5): 135-146. [本文引用:1]
[28] Schuster G T, Huang Y, Dai W, et al. Theory of multisource crosstalk reduction by phase-encoded statics[J]. Geophysical Journal International, 2011, 184(2): 1289-1303. [本文引用:1]
[29] Dai W, Fowler P, Schuster G T. Multi-source least-squares reverse time migration[J]. Geophysical Prospecting, 2012, 60(4): 681-695. [本文引用:1]
[30] Bai J, Chen G, Yingst D, et al. Attenuation compensation in viscoacoustic reverse time migration[C]//Expand ed Abstracts of the 83rd SEG Annual International Meeting, 2013: 3825-3830. [本文引用:1]
[31] Carcione J. Wave fields in real media: wave propagation in anisotropic, anelastic and porous media[M]. Elsevier Ltd: Pergamon Press, 2001. [本文引用:1]
[32] Lou M, Rial J A. Modelling elastic wawe propagation in inhomogeneous anisotyopic media by the pseudospectral method[J]. Geophysics, 120(10): 60-72. [本文引用:1]
[33] Carcione J M. A generalization of the Fourier pseudospectral method[J]. Geophysics, 2010, 75(6): A53-A56. [本文引用:1]
[34] Liu F Q, Zhang S, Morton A, et al. Anti-dispersion wave equation for modelling and reverse-time migration[C]//Expand ed Abstracts of the 78th SEG Annual International Meeting, 2008: 2277-2281. [本文引用:1]
[35] Lindeberg T. Scale-space theory in computer vision[M]. Kluwer Academic Publisher, 1994. [本文引用:1]
[36] Beydoun W B, Mendes M. Elastic ray-born l2-migration/inversion[J]. Geophysical Journal International, 1989, 97(1): 151-160. [本文引用:1]