2. 中国科学院大学, 北京 100049
2. University of Chinese Academy of Sciences, Beijing 100049, China
黏弹性叠前偏移方法通过在偏移过程中补偿地球介质对地震波的吸收,恢复被衰减的高频成份并校正频散,可获得较常规偏移方法更高分辨率的地下构造图像.由于这一过程遵循了地震波传播的物理规律,因此获得的高频成份是真实的,高分辨率成像结果也更能反应地下构造的实际情况.但在应用黏弹性叠前偏移方法进行高分辨率地震成像时,需利用Q值模型;由于Q值建模需利用地震信号随频率变化的幅值(Tonn, 1991),因此采用类似于深度域速度建模的方法进行Q场建模将遇到极大的困难.此外,当利用反射地震资料进行Q值建模时,除了获得准确的随频率变化的幅值十分困难外,另一个主要问题是地震反射的薄层调谐.薄层调谐是由于一组相近界面的反射波相互作用产生的效果,这在实际反射地震资料中是普遍存在的.薄层调谐将导致反射波的频谱发生较大的改变,这一改变甚至远大于吸收衰减的效果,因此不能简单地通过频谱改变决定实际介质的Q值(Zhang et. al,2013).
现行获得地球介质Q值的方法更多的是利用VSP数据和井间地震资料.由于VSP数据和井间地震资料中存在幅值占优的透射波,利用这些透射波信息可以获得透射波所穿过区域的Q值.但实际地震勘探中,并不总是进行VSP和井间观测;即使进行,数量也非常少,不足以描述勘探区域中Q值的横向变化.此外,尽管Q值是描述地震波吸收衰减的物理量,但在Q值建模时,还需考虑地震资料的信噪比;当高频噪音较发育时,若采用真实的Q值进行吸收衰减补偿,将放大高频噪音,极大地降低偏移成像的信噪比,反而不能实现提高分辨率的目标.因此,就黏弹性叠前偏移的Q值建模而言,利用VSP数据和井间地震资料也并不是一个理想的解决方案.
黏弹性叠前时间偏移(Zhang et. al, 2013, 2016)是新近提出的一种在偏移过程中补偿地球介质的吸收衰减的叠前偏移方法.不同于黏弹性叠前深度偏移方法(Mittet et al., 1995; Zhang and Wapenaar, 2002; Deng and McMechan, 2007; Wang, 2008; Xie et al., 2009; Zhu et al., 2014)使用每个成像点上的Q值,黏弹性叠前时间偏移使用一种新提出的等效Q值(Zhang et. al,2013).由于特定空间位置上的吸收补偿和频散校正仅由该位置处的等效Q值决定,而修改任一空间位置上的等效Q值,仅影响该位置处成像的吸收衰减补偿效果,可采用扫描等简单方法确定各样点处的等效Q值;这一点类似于叠前时间偏移中均方根速度建模.等效Q值的提出和应用,解决了黏弹性叠前偏移实际应用中面临的Q值建模的巨大困难,通过进一步引入仅对菲涅尔带内地震信号进行叠加成像的稳相偏移方法,有效压制粘性吸收补偿导致的噪音放大和偏移噪音,实现高分辨率、高信噪比叠前偏移成像(Zhang et al., 2016).这一方法已在不同区域的三维实际资料的应用中见到了很好效果.
现行的黏弹性叠前时间偏移利用等效Q值和均方根速度场.当成像点与炮、检点的水平距离较大且速度纵向变化较强时,基于均方根速度的走时计算存在误差,导致陡倾角构造和断层的断面不能很好成像.这一问题在常规的“直射线”叠前时间偏移方法中也存在.相对于常规偏移而言,黏性偏移的计算量比较大,通常是通过GPU加速计算,而GPU的固有性质是代数计算快而内存读写慢,这就使得常规的基于走时表的方法在GPU上不能得到较好的加速.本文针对黏弹性叠前时间偏移GPU算法实现的这一特点,通过在均方根速度基础上增加两个无量纲修正系数,改进了黏弹性叠前时间偏移的走时计算方法;这两个修正系数可由均方根速度场计算得到,也可根据成像结果进行修正.改进的走时计算有效提升了黏弹性叠前时间偏移对陡倾角构造和断层的成像效果.与基于走时表的“弯曲射线”叠前时间偏移方法相比,本文方法更适用于偏移方法的GPU加速(李博, 2009;刘国峰等,2009),避免了频繁读取走时表导致的计算时间增加(Wilt, 2013; John et al., 2014).
为补偿高频衰减和校正频散,黏弹性叠前时间偏移采用了频率域计算方法.频率域计算使得计算量增加了很多.为减少成像过程中频率域累加的计算量,本文提出了一种分时段的频率域算法.这一方法通过对记录的时域地震信号进行分段,增大频率间距,从而减少参与累加的频率域分量的数量.若将地震数据分为两段,理论上,就可将成像计算的计算量减少近一半.
本文将这一改进的黏弹性叠前时间偏移方法应用于实际三维陆上地震资料,并与商业处理系统的偏移成像效果进行对比.与商业系统的以及现行的黏弹性叠前时间偏移方法的成像结果对比表明,本文方法在获得高分辨率的成像结果的同时,明显改善了黏弹性叠前时间偏移方法对断层和陡倾角构造的成像效果,大幅提升了这一方法的计算效率.
1 黏弹性叠前时间偏移及等效Q值对任一地震道f(t),三维黏弹性叠前时间偏移的脉冲响应可表达为
(1) |
式中ω是角频率,ω0是地震道的主频,F(ω)是地震道f(t)的傅里叶变换,τs和τg是由均方根速度Vrms求得的炮点(xs, ys)和检波点(xg, yg)到成像点(x, y, T)的走时,
(2) |
Qeff就是等效Q值,与均方根速度类似,它是替代上覆地层各个不同Q值影响的一个等效参数,等效Q值和均方根速度可共同表达为
(3) |
式(3)中Ql和vl为上覆各层介质的Q值和速度,ΔTl是各层介质的垂直旅行时厚度.
式(1)的单地震道偏移成像公式表明,任一成像点的黏弹性偏移成像,仅与该成像点处的均方根速度Vrms和等效Q值Qeff有关,而修改任一空间位置上的均方根速度或等效Q值,仅影响该位置处成像的聚焦和吸收衰减补偿效果,因此,可以应用扫描的方法确定任一空间位置上的均方根速度和等效Q值.式(3)表明,黏弹性叠前时间偏移使用了与常规叠前时间偏移相同的均方根速度,因此该方法可直接应用常规叠前时间偏移的偏移速度(这也是命名为黏弹性叠前时间偏移的原因).而所谓扫描确定等效Q值,就是令该空间位置上的等效Q值取为一系列可能的值,对比不同数值的吸收衰减补偿效果(高频恢复情况),最终确定一个最佳的等效Q值.等效Q值的引入,极大地简化了Q值建模,克服了黏弹性叠前深度偏移面临的Q值建模的巨大困难.补偿吸收衰减偏移方法会带来高频端的噪音问题,通过给定补偿阈值和补偿高截止频率来控制高频端的不稳定性.
基于倾角道集,获得在各成像点用时间域倾角表达的二维菲涅尔带(φx±, φy±), 可得三维黏弹性叠前时间偏移的成像公式为( Zhang et.al,2016)
(4) |
式中下标m是地震道的序号,H是地震道的偏移距,循环叠加是对全部地震道按偏移距大小进行, n是地震道总数;Ω(x, y, T, φx±, φy±, θx, θy)是各成像点上由二维菲涅尔带决定的权系数,其中(θx, θy)是每个地震道的入射、反射射线在成像点处所对应的(虚拟)反射面的倾角,可由下式计算得到
(5) |
若满足φx-<θx<φx+和φy-<θy<φy+, Ω=1,否则Ω=0.权系数Ω(x, y, T, φx±, φy±, θx, θy)使式(4)的黏弹性叠前时间偏移算法实现了时-空变的最优偏移孔径(仅对菲涅尔带内地震信号进行叠加成像),避免了常规偏移孔径不能兼顾多个构造特征的问题,其在倾角域反假频的同时,能够直观的展示对成像贡献的菲涅尔带,基于菲涅尔带的叠加,可以有效压制粘性吸收补偿导致的噪声放大和偏移噪声.
借用常规叠前时间偏移的偏移速度场,根据Q值扫描确定等效Q值,利用倾角道集估计各成像点的菲涅尔带,就可用式(4)的偏移公式进行黏弹性叠前时间偏移成像,减少了菲涅尔带之外噪声对成像结果的影响,获得高分辨率、高信噪比的CRP成像道集.实际应用中,整个三维成像空间上的等效Q值和菲涅尔带可利用选定样点处的结果进行插值得到.
2 高精度走时计算本节讨论如何基于黏弹性叠前时间偏移的GPU算法改进走时计算方法,以提高该方法对断层和陡倾角构造的成像精度.式(2)所式的走时计算公式是利用如下近似式得到的(Zhang et al., 2016)
(6) |
当式(6)中介质速度vl沿时间深度方向变化较剧烈,即
解决这一问题的一个共同方法是在各个横向位置(CDP)处,假设介质为水平层状,利用射线法等方法计算炮、检点到成像点的走时;而该横向位置处水平成层介质的层速度,是由均方根速度反演(如Dix公式)得到.就常规叠前时间偏移而言,其实现方法更多采用的是走时表算法,即事先对每一条成像线,计算好由CDP坐标、炮(检)点到成像点的水平距离、时间深度这3个维数决定的走时表,偏移计算时由此表插值确定炮(检)点到成像点的走时.需指出的是,这类基于水平层状介质假设的走时算法是不能直接用来确定(或修正)偏移速度的,它需在利用式(2)的计算公式,由扫描等方法确定了光滑的均方根速度后,才能应用于偏移计算.
由于黏弹性叠前时间偏移必须通过频率域积分实现,计算量相对常规成像而言要大很多,常规CPU计算能力很难满足实际生产需要,必须使用GPU来加速偏移,而走时表方法涉及了巨大走时表的频繁读取,因此不适应GPU加速计算.为此,本文提出了一个包含(炮点或检波点到成像点的)水平距离高次项的高精度走时算法.不同于包含水平距离高次项的多项式近似(Taner and Fulton, 1969),本文通过在均方根速度基础上增加两个无量纲修正系数,改进了式(2)的走时计算,具体如下:
(7) |
式中
当α=0和β=0,式(7)即退化为式(2),无量纲参数α和β与均方根速度Vrms同样,仅是成像点的(x, y, T)的函数,改变任一空间位置上的α和β,仅影响该位置处大入(出)射角度的偏移同相轴的聚焦.
虽然我们提出在黏弹性叠前时间偏移GPU算法实现中以增加两个无量纲参数α和β来替代走时表策略改进走时的计算精度,但对于任一成像点的无量纲参数α和β的求取,仍可借鉴走时表方法,即基于均方根速度反演得到的对应该成像点横向位置的水平成层速度模型,用最小二乘方法求得.具体方法是:在水平成层的二维速度模型中,从成像点按等间距变化的出射角度,向地表发出一系列射线(见图 2),计算地表接收点的水平距离和走时,将一系列水平距离和走时的数据对代入式(7),可由最小二乘拟合求得该成像点的α和β.在(7)式中将α和β作为优化系数可以得到(8)式(由于(7)炮点走时和检波点走时公式类似,此处只给出炮点公式):
(8) |
将(8)式写成矩阵形式得到(9)式:
(9) |
对于不同角度出射构成的超定方程,利用超定方程的最小二乘公式ATAm = ATb可以求得α和β.等间距变化的出射角度使得小距离的数据对多,而大距离的数据对少,这相当于给式(7)的拟合引入了一个随水平距离减少的权系数,从图 2中等角度间隔出射的射线可以看到在到达地表后,间隔会随着角度的增大而变大.该权重可以很好地保持近距离端的稳定性,又具有大距离端的修正能力,使得陡倾角构造较好成像,又不影响其他构造的成像.
除利用水平层状速度模型求取无量纲参数外,也可通过扫描方法,由陡倾角构造的最佳成像,确定和修改无量纲参数α和β.这弥补了走时表方法不能调整偏移速度的缺点,从而可提升本文黏弹性叠前时间偏移方法对陡倾角构造的成像效果.
图 1给出了实际地质模型在一个典型CDP点上的均方根速度随时间深度变化的曲线和基于该曲线由Dix公式反演得到的时间域水平成层速度模型,图中Vrms表示该CDP上的均方根速度,Vint指示该CDP上反演得到的时间域层速度.图 2给出了图 1模型中成像点A向地表发出的系列射线的示意图.就成像点A、B和C,利用最小二乘拟合方法,求得的无量纲参数α和β为(0.0213,0.004)、(0.0345,-0.0028)和(0.0358,-0.004). 图 3显示了成像点A、B和C处,式(7)和式(2)计算的走时与基于射线追综的理论走时对比的走时误差.从图中可看出,常规走时在入(出)射角度较大时和基于射线追踪的走时差异较大,在三个点处分别为12 ms, 43 ms和54 ms.通过引入两个无量纲参数,式(7)的走时计算公式明显改善了大入(出)射角度方向传播波场走时计算的精度,在70°角度范围内保证了走时误差小于0.07 ms.这意味着当地下介质横向速度变化较平缓时,应用式(7)可对倾角达到70°的构造正确成像.
黏弹性叠前时间偏移可获得较常规叠前时间偏移更高分辨率的成像剖面.然而,与有着很高计算效率的常规叠前时间偏移相比,这一方法实现环节包含的频率域积分产生了巨大的计算量,除上文提到的需基于GPU加速实现其算法外,本节同时提出一种分时段的频率成像策略以进一步提高计算效率.
当1/Qeff=0,式(1)可简化为
(10) |
式中f′(t)是地震道f(t)的一阶时间导数.式(8)就是常规叠前时间偏移的脉冲响应,对比式(10)和式(1)可知,相对于有着很高计算效率的叠前时间偏移方法,黏弹性叠前时间偏移增加的计算量来自于式(1)的频率域积分,但这一频率域积分(或采用时间域褶积)是补偿黏性吸收和校正频散必需的.为提高黏弹性叠前时间偏移的计算效率,引入式(5)的倾角约束,通过减少需要进行式(1)或(4)的频率域积分的成像点数量,继而减少了需要累加的面积,以此来降低偏移的计算量.
实际偏移时,式(1)或(4)的频率域积分是采用下式的频率域累加实现的:
(11) |
式中Δω=2π/T0,T0是地震道地震记录的时长;l1=ωmin/Δω,l2=ωmax/Δω,ωmin和ωmax是地震信号的有效低频和高频,函数
将地震记录的时长等分为两段,每段的时长是T0/2+TU和T0/2+TD,其中TU和TD是两段的重叠部分,则Δω将增加近一倍,其取值将在下面讨论.在偏移计算过程中,当τs+τg≤T0/2时,将前半时段地震记录的傅里叶变换代入式(11)进行计算.则仅需近似一半的计算量就可完成式(11)的偏移计算;当τs+τg>T0/2时,将后半时段地震记录的傅里叶变换和新的τs+τg=τs+τg-T0/2+TD代入式(11),同样仅需近似一半的计算量就可完成式(11)的偏移计算.若地震记录的时长很长,则可将地震记录等分为三段,此时计算时间将近似为原来的三分之一.由于重叠部分的时长是由(τs+τg)/Qeff和函数的最大值决定的所分段数n并不是越大越好,它是由T0/n与TU和TD的大小对比决定的.若两者相等,分段计算将不减少计算量,且增加了计算的复杂度.下面将讨论TU和TD的取值.
利用傅里叶变换的褶积定理,式(1)改写为时间域褶积:
(12) |
式中函数h是一个由参数c=(τs+τg)/Qeff和函数ϕ的最大值决定的时变子波,有
(13) |
若地震记录的时长等分为两段,则将τs+τg=T0/2和对应空间位置的最小Qeff代入式(13),可求得相应的子波hB,根据子波hB的宽度,可求得对应的TU和TD.当1/Qeff=0,函数h退化为脉冲函数,式(13)变得和式(10)一样.
图 4给出了T0=4 s时将地震记录等分为两段的子波hB,其中τs+τg=2 s,选取的Qeff=90,函数ϕ是最大值为100的e指数函数(Zhang et al., 2013).在图 4中,将小于最大振幅1%部分截断,局部放大结果见图 5,其中虚线表示最大振幅1%的界限,可以求得TU=168 ms,TD=110 ms.图 6给出了一个实际地震道采用式(9)得到的一个偏移成像道和采用分两段的分段算法得到的偏移成像道的局部对比,图中显示了两种算法的计算结果.从图 6可知,采用分段算法保持了计算的精度.
为了验证本文方法的有效性,本文给出我国东北部的三维实际资料的应用效果.该工区采用14L×28S×168R-45°斜交观测,面元大小为25 m×25 m,覆盖次数为12(纵)×7(横)=84次,接收道为2352,道距为50 m.该区块目的层为2 s到3 s,该工区构造复杂、断层多,而且断块下的同相轴接触关系复杂,提高分辨率和构造精准成像对其是关键.图 7是改进走时后的某条成像线的成像结果.图 8a为常规商业软件偏移结果,图 8c为常规走时的高分辨率结果.对比3幅图可以看到,改进后的高分辨率成像结果在构造上和商业软件结果一致.其分辨率比商业软件高,细节更清晰.改进后的高分辨率成像结果和改进前相比可以看出,其对陡倾角构造成像更好,对断层断块的刻画更为清晰.图 9进一步给出了图 8a和8b的叠加频谱对比,实线是商业软件偏移结果,虚线是本文方法偏移结果.从图 9中可以看出,采用相同的输入数据,不需要额外任何资料,本文的黏弹性偏移成像结果较商业软件处理结果有效频带拓宽了16 Hz.
针对图 7的一条测线的黏弹性偏移成像(满覆盖叠前资料250G),分别应用了不分段(常规方法)和本文的分两段计算方法,在完全相同的硬件环境下(GTX980Ti),两者的计算耗时分别是1870 s和1107 s,见表 1.考虑到读取数据的耗时,采用分两段黏弹性偏移计算方法,可整体提高计算效率1/3.
本文修正了适合于GPU加速的黏弹性叠前时间偏移走时计算方法,改进了该方法对陡倾角构造和断层的成像效果.这一计算方法,针对GPU代数运算快的特点,运用多参数描述方法来取代常规商业软件中的走时表方法,充分发挥了GPU计算能力,规避了频繁内存读取引发的低计算效率问题.本文还进一步提出了分时段的频率域成像策略,在不影响成像精度的情况下,减少了黏弹性叠前时间偏移固有的频率域积分计算强度,大幅提升了黏弹性叠前时间偏移的计算效率.本文将改进的黏弹性叠前时间偏移方法应用于三维陆上地震资料.与现行的商业偏移软件对比表明,该方法不仅获得了更高分辨率的成像结果,也实现了对断层和陡倾角构造的清晰成像,而新方法的计算耗时也较改进前减少了三分之一以上,对大规模工业应用有很重要的意义.
Deng F, McMechan G A.
2007. True-amplitude prestack depth migration. Geophysics, 72(3): S155-S166.
DOI:10.1190/1.2714334 |
|
John C, Grossman M, McKercher T. 2014. Professional CUDA C Programming. Wrox Pr/Peer Information Inc.
|
|
Li B, Liu G F, Liu H.
2009. A method of using GPU to accelerate seismic pre-stack time migration. Chinese Journal of Geophysics, 52(1): 245-252.
|
|
Liu G F, Liu H, Li B, et al.
2009. Method of prestack time migration of seismic data of mountainous regions and its GPU implementation. Chinese Journal of Geophysics, 52(12): 3101-3108.
|
|
Mittet R, Sollie R, Hokstak K.
1995. Prestack depth migration with compensation for absorption and dispersion. Geophysics, 60(5): 1485-1494.
DOI:10.1190/1.1443882 |
|
Taner M T, Fulton K.
1969. Velocity spectra-digital computer derivation and applications of velocity functions. Geophysics, 34(6): 859-881.
DOI:10.1190/1.1440058 |
|
Tonn R.
1991. The determination of the seismic quality factor Q from VSP data:A comparison of different computational methods. Geophysical Prospecting, 39(1): 1-27.
DOI:10.1111/gpr.1991.39.issue-1 |
|
Wang Y H.
2008. Inverse-Q filtered migration. Geophysics, 73(1): S1-S6.
DOI:10.1190/1.2806924 |
|
Wilt N. 2013. The CUDA Handbook: A Comprehensive Guide to GPU Programming. Amsterdam: Addison-Wesley Professional.
|
|
Xie Y, Xin K F, Sun J, et al. 2009. 3D prestack depth migration with compensation for frequency dependent absorption and dispersion. //79th Annual International Meeting, SEG. Expanded Abstracts, 2919-2922.
|
|
Zhang J F, Wapenaar K.
2002. Wavefield extrapolation and prestack depth migration in anelastic inhomogeneous media. Geophysical Prospecting, 50(6): 629-643.
DOI:10.1046/j.1365-2478.2002.00342.x |
|
Zhang J F, Wu J Z, Li X Y.
2013. Compensation for absorption and dispersion in prestack migration:An effective Q approach. Geophysics, 78(1): S1-S14.
DOI:10.1190/geo2012-0128.1 |
|
Zhang J F, Li Z W, Liu L N, et al.
2016. High-resolution imaging:An approach by incorporating stationary-phase implementation into deabsorption prestack time migration. Geophysics, 81(5): S317-S331.
DOI:10.1190/geo2015-0543.1 |
|
Zhu T Y, Harris J M, Biondi B.
2014. Q-compensated reverse-time migration. Geophysics, 79(3): S77-S87.
DOI:10.1190/geo2013-0344.1 |
|
李博, 刘国峰, 刘洪.
2009. 地震叠前时间偏移的一种图形处理器提速实现方法. 地球物理学报, 52(1): 245–252.
|
|
刘国峰, 刘洪, 李博, 等.
2009. 山地地震资料叠前时间偏移方法及其GPU实现. 地球物理学报, 52(12): 3101–3108.
DOI:10.3969/j.issn.0001-5733.2009.12.019 |
|