2011, Vol. 54
2. 中国科学院地质与地球物理研究所, 北京 100029
2. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
为了适应深层隐蔽油气藏勘探的需求，特别是提高对薄互层、小砂体、复杂小断层等深层复杂地质体的识别能力，迫切需要发展一种高分辨率叠前地震成像技术，使其能够在成像中直接补偿地震波的高频能量损失[1~8]，以获取高分辨率反射系数成像，为我国陆相深层隐蔽油气藏勘探奠定坚实的技术基础[9].
工业界通常采用零相位反褶积和时变谱白化、反Q滤波[10~16]、黏性介质叠前深度偏移[17~19]等方法补偿介质的黏性吸收作用.零相位反褶积和时变谱白化虽然可补偿地震波的频率损失，却无法校正地震波的时变相位扭曲；反Q滤波方法虽然可克服前二者的缺点，但叠后反Q滤波无法正确弥补叠加时的高频损失，叠前反Q滤波又无法处理非均匀Q值；黏性介质叠前深度偏移方法可处理横向剧烈变化的Q值场，但巨大的计算量以及在非均匀Q值建模方面的困难，制约了其在实际数据处理中的应用.
频域黏性介质叠前时间偏移技术[20]通过定义等效的叠加Q因子和在偏移过程中扫描确定叠加Q值，实现了直接在叠前时间偏移成像中补偿地震波黏性吸收作用.该方法通过控制最大的补偿因子，保证了黏性吸收补偿的稳定性.通过保证浅、中、深层具有相同的频带，求得了等效的叠加Q值.克服了反Q滤波和黏性叠前深度偏移在叠前资料处理中的缺陷，在对中等复杂速度构造实现准确偏移成像的同时补偿黏性吸收，有效地提高了深层成像的分辨率.
频域算法涉及五维数据体的计算，计算量仍十分巨大，没有充分体现出叠前时间偏移的效率优势.因此，进一步提高黏性介质叠前时间偏移的计算效率是实现海量数据处理的关键.显式空间褶积短算子是叠前深度偏移中为了提高计算效率而广泛使用的一种设计方案[17~19, 21~26]，本文将这一优化设计思想引入到黏性介质叠前时间偏移中，采用加权最小平方优化设计方案获取时域最优褶积短算子，将频域补偿方案转化为时域补偿方案，将频域大量的频率相乘运算转化为少量的时域褶积运算，减少一个处理循环维，借以提高黏性介质叠前时间偏移的计算效率，而这方面的研究工作尚无文献涉及.
2 时域算法的基本原理在二维情况下，由单程波方程和稳相点原理出发，推导出了频率-空间域黏性介质相移补偿公式的高频渐进解为[20, 27~29]
|
(1) |
其中:j为虚数单位；ω 为角频率；T为总的垂直走时；x为炮(检波)点到成像点的水平距离；crms为均方根速度；P(x，ω，T)为频率-空间域的地震波场；f(ω)为输入道的频谱；px0为反射点所对应的射线参数；Γ(ω，px0)为补偿因子项，其表达式为
|
(2) |
其中，cQrms 为与Q有关的均方根速度；Qeff为叠加品质因子；设A为公式(1)中的幅值项，则
|
(3) |
对于给定任一道原始地震数据，则可知该数据道与震源之间的距离xs 和检波点之间的距离xr；对任一拟成像的CDP点和成像的时间深度T，根据公式(1)计算出震源正传和单道数据反传的幅值和走时，分别记为As、Ar 和ts、tr.利用波动方程深度偏移的反褶积成像条件对输入道数据成像，若设震源为一时间脉冲，有
|
(4) |
其中，I(x，T)为黏性介质叠前时间偏移的成像；在积分中的前三项可看作是该道数据求取半导数后的傅氏变换，其在时间域记为Dh(t)，即
|
(5) |
其中:F-1 代表反Fourier变换.根据Fourier变换性质，频域相乘对应于时域的褶积，将式(4)变换到时间域可进一步简化为
|
(6) |
其中，t代表地震波的走时，且t=ts+tr；S(t)为时域长褶积算子，其Fourier变换等于频域补偿因子谱.式(6)表明在时间域对原始数据求取半导数，然后与时域褶积算子S(t)在ts+tr 时刻褶积成像，再乘上由幅值确定的权系数，即可获取最终的成像结果，这是时域黏性介质叠前时间偏移的基本原理.由于数据是离散的，上述幅值的拾取是通过两点插值实现的.
令频域补偿因子谱为W′(ω)，则
|
(7) |
|
(8) |
由于S(t)很长，为了提高计算效率，必须采取一种优化策略，将S(t)截短，并使其频谱在有效频带内能够精确拟合补偿因子谱，从而避免对S(t)直接截短所产生的吉布斯现象.
3 稳定性控制由于补偿因子在补偿地震波的高频能量损失的同时，也会增强数据中的数值误差或高频噪声而使算法表现出固有的不稳定性[30].因此，时域褶积短算子算法的稳定性是时域黏性介质叠前时间偏移中最为关键的问题.在短算子设计过程中，稳定性控制必须同时满足压制高频噪声和提高短算子设计精度这两项基本要求.因此，本文提出了一种保持高频的、门限阈值约束的、时变的稳定性控制算法.
该稳定性控制方法的实质是在频域内通过控制最大补偿因子W′(ω)，将其改造成稳定的补偿因子谱W(χ)(χ =ωt，χ 与走时t有关，因而W(χ)是时变的)，即:设定时变的频率截止值χc(χc 值由幅度控制门限阈值Glim计算得到).当χ＜χc 时，W(χ)为精确的反Q补偿因子谱；当χ≥χc 时，进行高频压制，利用一个一阶光滑可导的函数α(χ)改造补偿因子谱W′(ω)，使其光滑地衰减到给定的阈值Glim, 即
|
(9) |
该方法的优点是可以通过改变阈值大小来控制时变补偿频带范围，使其既可以保持稳定的补偿效果又尽可能地保持有效的高频信息.
4 时域褶积短算子设计采用加权最小平方方法设计时变褶积短算子的核心思想是:在时间域内设计一个尽可能短的一维时域补偿褶积短算子U(t′)，使其傅氏变换后的频谱在有效频带范围内很好地拟合经稳定性控制改造后的反Q补偿因子谱W(ω).这个问题可以由 Fourier变换表达成如下形式:
|
(10) |
其中:t′为时间变量，将式(10)写成离散形式:
|
(11) |
其中:M为W(ω)在频域内的采样点数，N为设计短算子的长度.为了提高计算效率，N一般应少于30个采样点比较适宜.
对于零相位子波，其关于t′=0对称，则短算子的频域响应可以通过对时域短算子做余弦变换得到，即
|
(12) |
令C为余弦变换的系数矩阵，U代表短算子向量，W代表频域反Q补偿因子向量，则有
|
(13) |
设加权预测误差函数用ε 表示:
|
(14) |
其中:$\tilde{e}$=C$\tilde{U}$-W，$\tilde{U}$ 为U的一个最小平方近似.H代表复共轭转置.Λ 为权重系数的对角阵.权系数给出准则是:在有效频带范围内给予较高的权重(设为1)，有效频带以外给予较低的权重(一般取为0.0005)[31].优化的目标是使预测误差函数ε 达到最小.令
|
(15) |
则推导得:
|
(16) |
该实数方程组可以采用共轭梯度或LU 分解等方法求解出最优褶积短算子.图 1、图 2分别给出了时域最优褶积短算子在频域和时域的表现形式；时域内，最优褶积短算子与稳定补偿因子谱经反Fourier变换到时域的理论长子波按短算子长度截取后的波形具有很好的一致性，说明构造的褶积短算子具有很高的精度；频域内，最优褶积短算子的频谱在有效频带内精确拟合了稳定补偿因子谱W(ω)，这也是时域算法精确补偿成像的根本保证；在有效频带外，短算子的频谱光滑地衰减下来，保证了算法的稳定性.因此，时域黏性介质叠前时间偏移就是通过该高精度的时变最优褶积短算子$\tilde{U}$ 与输入数据道相褶积来完成的.
|
图 1 时域短算子与长子波的对比图中实线代表将稳定补偿因子谱反变换到时域后长子波截取波形，虚线代表时域短算子. Fig. 1 The comparison of short operator in time domain with long wavelet The solid line represent for cutoff short wavelet transformed from spectrum after smoothing in frequency domain.The dash line represent for short operator calculated by weighted least squares method. |
|
图 2 短算子的频谱与稳定补偿因子谱的对比 图中实线代表稳定补偿因子谱，虚线代表短算子经Fourier变换后的频谱. Fig. 2 The comparison of short operator spectrum with spectrum after smoothing The solid line represent for the spectrum after smoothing in frequency domain.The dash line represent for spectrum of short operator. |
由公式(7)可知，褶积短算子的补偿量只与走时t和Qeff比值有关，在实践中采用表驱动策略，事先按照t/Qeff 等间距变化建好短算子系数表库，将短算子设计与时域补偿偏移运算相剥离；计算时可按照t/Qeff 依据表的间距取整，到表中直接调取短算子系数即可参与运算.短算子系数表一旦建好，可以重复使用，大大地提高了算法的计算效率.
5 偏移流程对经过叠前去噪和静校正后的叠前数据进行时域黏性介质叠前时间偏移，需要以下几个步骤:
(1) 根据待处理数据共深度点(CDP)的分布范围确定成像空间的大小，给出最大偏移孔径成像角度及衰减系数等必要的成像控制参数.对叠前地震数据求取半导数；计算走时、振幅及短算子系数表.
(2) 选定几个共中心点(CMP)道集作动校正，初步确定一个初始的横向均匀叠加速度场c0.
(3) 利用c0 对原始数据进行常规叠前时间偏移，得到共反射点(CRP)道集，再用c0 对其反动校，然后再作动校正将CRP道集校平，将该动校正速度作为该点新的叠加速度.
(4) 对新得到的叠加速度场平滑处理后作为偏移的速度场c1.利用c1 对原始数据进行常规叠前时间偏移，得到校平后CRP 道集，叠加后便得到准确构造成像剖面.
(5) 采用改进对数谱比法扫描分析叠加Q值(详细方法另文论述)，获取准确的叠加Q值场，对其平滑后作为黏性偏移的输入.
(6) 采用基于短算子的时域黏性介质叠前时间偏移方法，对于每一个成像点，计算出该点的t/Qeff, 查表得到短算子系数，与求取半导数后的地震数据按照走时ts+tr 褶积成像，乘上相应的权系数，得到补偿校平后的CRP道集，叠加后获取高分辨率的构造成像剖面.补偿过程中增强的随机噪声会在叠加中得到很好的压制，这也是叠前反Q补偿的主要优势之一.
6 数值算例首先，考察时域和频域算法在补偿效果的一致性及效率改善情况.采用相移法正演模拟了具有三个水平反射界面的均匀黏性介质的炮记录数据(界面处人为设定了反射系数)，其中，介质模型速度为2000m/s, Q值为50；反射界面之间的厚度均为1000m.震源子波的主频为30Hz, 记录的采样点数为2000，采样间距为0.002s, 道间距为10 m, 单炮数据中共包含500道数据.
从正演模拟的炮数据中抽取一道近偏移距数据，分别采用频域算法和时域算法获取该道的偏移脉冲响应(图 3)，对比发现二者的补偿效果几乎一致；然后再利用这两种算法对该炮记录进行偏移成像(图 4)，两个偏移剖面也取得了几乎一致的补偿效果，深层反射界面的分辨率与浅层相一致.在计算效率方面，常规叠前时间偏移的处理时间为116s；频域算法所用的时间为1085s；时域算法(短算子长度为25)所用时间为138s, 其计算效率比频域算法提高了近8倍，与常规叠前时间偏移的计算效率大体相当，这说明基于最优褶积短算子的时域算法可极大地提高计算效率.
|
图 3 (a)采用频域算法得到的偏移脉冲响应；(b)采用时域算法得到的偏移脉冲响应 Fig. 3 Migration impulse responses for (a) frequency domain algorithm and for (b) time domain algorithm |
|
图 4 水平层状理论模型偏移结果 (a)频域算法；(b)时域算法. Fig. 4 The migration profile of horizontally stratified model (a) The algorithm in frequency domain； (b) The algorithm in time domain |
其次，考察时域算法在复杂地质模型情况下的适应能力.采用二维频域有限差分方法利用双程声波波动方程建立黏性介质的衰减正演模拟理论模型.正演模拟所用的宏观介质参数模型如图 5所示.模拟后的数据包含51个共炮点道集，震源采用的是最大频率为55Hz的雷克子波，炮点分布位置是从447.75~3147.75 m, 炮间距54 m；每炮800 道，检波点间距为4.5m；记录长度1.68s, 采样间隔4ms；在模型两侧每边共设20 个网格点的衰减区以减少边界反射.为了保证模拟数据更接近真实数据，在模拟中保留了表面多次波和层间多次波.
|
图 5 用于正演模拟所用的宏观介质参数模型 Fig. 5 Subsurface real phase velocity and Q model |
分别采用频域算法和时域算法对上述理论模型数据进行补偿偏移，采用共反射点叠加方案，获取补偿成像剖面(图 6).从图 6上看，所有的反射轴均正确归位并精确成像，断面清晰；右侧楔状体的上界面处，地震波从高速介质入射到低速介质，波形出现极性反转.这说明无论是频域算法还是时域算法都可以对复杂构造精确成像.此外，在对黏性吸收给出正确补偿后，深层和远偏移距的子波频带得到拓宽，空间子波被压缩，振幅增强，深层获得了与浅层相一致的分辨率.该数值算例说明，本文给出的时域算法可以改善地震成像剖面的分辨率，且时域算法拥有与频域算法相同的补偿效果.
|
图 6 (a)频域的黏性补偿叠前时间偏移剖面；(b)基于时域褶积短算子的黏性补偿叠前时间偏移剖面 Fig. 6 The prestack time migration profile of compensation for anelastic absorption (a) in frequency domain and (b) basing on short convolution operator in time domain |
(1) 本文给出了基于叠加Q值的时域黏性介质叠前时间偏移算法，其补偿量只与成像点处的叠加Q值有关，因而允许每一道数据在同一时刻的Qeff不同，即Qeff可在深度方向连续变化，横向上弱变化，这比频域层状反Q滤波算法具有更广的适用性，可应用于叠前地震数据的高分辨率成像.
(2) 本文通过控制最大补偿因子，保证了黏性吸收补偿的稳定性.采用加权最小平方优化计算方案获取时域最优褶积短算子，将大量的频率相乘运算转化为少量的时域褶积运算.同时采取表驱动策略，事先建好短算子系数表库，将褶积短算子的计算过程与补偿成像过程相剥离，大大提高了算法的计算效率.因此，基于短算子设计的时域黏性介质叠前时间偏移可以很好地结合到现有的常规叠前时间偏移方案中.
(3) 理论正演模拟数据处理结果表明，本文给出的算法可以改善地震成像剖面的分辨率，时域算法拥有与频域算法相同的补偿效果.但计算效率得到了极大地改善，与常规叠前时间偏移的计算效率相接近.
| [1] | Robinson J C. A technique for the continuous representation of dispersion in seismic data. Geophysics , 1979, 44(8): 1345-1351. DOI:10.1190/1.1441011 |
| [2] | Zhao B, Zhou H W, Li X G, et al. Attenuation analysis on commercial and low-saturation gas reservoirs. Antonio Ann. Internat Mtg., Soc. Expl. Geophys.. Expanded Abstracts , 2007: 1412-1416. |
| [3] | Causse E, Ursin B. Asymptotic error analysis of constant-velocity viscoacoustic migration. Geophysics , 1999, 64(4): 1036-1045. DOI:10.1190/1.1487043 |
| [4] | Kjartansson E. Constant Q wave propagation and attenuation. Journal of Geophysics Research , 1979, 82(B9): 4737-4748. |
| [5] | Carcione J M, Kosloff D, Kosloff R. Wave propagation simulation in a linear viscoelastic medium. Geophysics Journal International , 1999, 95(3): 597-611. |
| [6] | Emmerich H, Korm M. Incorporation of attenuation into time-domain computations of seismic wavefields. Geophysics , 1987, 52(9): 1252-1264. DOI:10.1190/1.1442386 |
| [7] | Grimbergen J L T, Dessing F J, Wapenaar C P A. Modal expansion of one-way operators in laterally media. Geophysics , 1998, 63(3): 995-1005. DOI:10.1190/1.1444410 |
| [8] | Keers H, Vasco D W, Johnson L R. Viscoacoustic cross-well imaging using asymptotic waveforms. Geophysics , 2001, 66(3): 861-870. DOI:10.1190/1.1444975 |
| [9] | Gray S H, Etgen J, Dellinger J, et al. Seismic migration problems and solutions. Geophysics , 2001, 66(5): 1622-1640. DOI:10.1190/1.1487107 |
| [10] | Bickel S H, Natarajan R R. Plane-wave Q deconvolution. Geophysics , 1985, 50(9): 1426-1439. DOI:10.1190/1.1442011 |
| [11] | Bickel S H. Constant Q inverse filters. 51st Ann. Internat. Mtg., Soc. Expl. Geophys.. Expanded Abstracts , 1981: 172-181. |
| [12] | Hargreaves N D, Calvert A J. Inverse Q filtering by Fourier transform. Geophysics , 1991, 56(1): 519-527. |
| [13] | Hargreaves N D. Similarity and the inverse Q filter: some simple algorithms for inverse Q filtering. Geophysics , 1992, 57(2): 944-947. |
| [14] | 姚振兴, 高星, 李维新. 用于深度域地震剖面衰减与频散补偿的反Q滤波方法. 地球物理学报 , 2003, 46(2): 229–233. Yao Z X, Gao X, Li W X. The forward Q method for compensation attenuation and frequency dispersion used in the seismic profile of depth domain. Chinese J. Geophys. (in Chinese) , 2003, 46(2): 229-233. |
| [15] | Wang Y H. A stable and efficient approach of inverse Q filtering. Geophysics , 2002, 67(2): 657-663. DOI:10.1190/1.1468627 |
| [16] | Wang Y H. Quantifying the effectiveness of stabilized inverse Q filtering. Geophysics , 2003, 68(1): 337-345. DOI:10.1190/1.1543219 |
| [17] | Mittet R, sollie R, Hokstad K. Prestack depth migration with compensation for absorption and dispersion. Geophysics , 1995, 60(5): 1485-1494. DOI:10.1190/1.1443882 |
| [18] | Mittet R. A simple design procedure for depth extrapolation operators that compensate for absorption and dispersion. Geophysics , 2007, 72(2): 105-102. |
| [19] | Zhang J F, Wapenaar K. Wavefield extrapolation and prestack depth migration in anelastic inhomogeneous media. Geophysical Prospecting , 2002, 50(6): 629-643. DOI:10.1046/j.1365-2478.2002.00342.x |
| [20] | 李雪英. 粘性吸收补偿:从反Q滤波到叠前时间偏移. 北京: 中国科学院研究生院, 2009 . Li X Y. Compensation for anelastic absorption: from inverse Q filtering to prestack time migration (in Chinese). Beijing: Graduate University of Chinese Academy of Sciences, 2009 . |
| [21] | Holberg O. Towards optimum one-way wave propagation. Geophysical Prospecting , 1988, 36(2): 99-114. DOI:10.1111/gpr.1988.36.issue-2 |
| [22] | Blacquière G, Debeye H W J, Wapenaar C P A, et al. 3D table-driven migration. Geophysical Prospecting , 2006, 37(8): 925-958. |
| [23] | Hale D. 3-D depth migration via McClellan transformations. Geophysics , 1991, 56(11): 1778-1785. DOI:10.1190/1.1442990 |
| [24] | Hale D. Stable explicit depth extrapolation of seismic wavefields. Geophysics , 1991, 56(11): 1770-1777. DOI:10.1190/1.1442989 |
| [25] | Maeland E. On the construction of the 3D band-limited extrapolation operator in the space-frequency domain. Geophysical Prospecting , 1993, 41(5): 645-658. DOI:10.1111/gpr.1993.41.issue-5 |
| [26] | Zhang J F, Verschuur D J, Wapenaar C P A. Depth migration of shot records in heterogeneous, transversely isotropic media using optimum explicit operators. Geophysical Prospecting , 2001, 49(3): 287-299. DOI:10.1046/j.1365-2478.2001.00255.x |
| [27] | 李雪英, 吕喜滨, 张江杰, 等. 稳定高效的时域反Q滤波方法. 地球物理学进展 , 2010, 25(1): 211–218. Li X Y, Lü X B, Zhang J J, et al. A stable and efficient approach of inverse Q filtering in time domain. Progress in Geophysics (in Chinese) , 2010, 25(1): 211-218. |
| [28] | 董春晖, 张剑锋. 起伏地表下的直接叠前时间偏移. 地球物理学报 , 2009, 52(1): 239–244. Dong C H, Zhang J F. Prestack time migration including surface topography. Chinese J. Geophys. (in Chinese) , 2009, 52(1): 239-244. |
| [29] | 董春晖. 叠前时间偏移:基于波动方程的架构、发展和应用. 北京: 中国科学院研究生院, 2008 . Dong C H. Prestack time migration: a wave equation based approach, development and practical aspects (in Chinese). Beijing: Graduate University of Chinese Academy of Sciences, 2008 . |
| [30] | Wang Y H. Inverse Q-filter for seismic resolution enhancement. Geophysics , 2006, 71(3): 51-60. DOI:10.1190/1.2192912 |
| [31] | Thorbecke J W, Wapenaar K, Swinnen G. Design of one-way wavefield extrapolation operators, using smooth functions in WLSQ optimization. Geophysics , 2004, 69(4): 1037-1045. DOI:10.1190/1.1778246 |