地球物理学报  2014, Vol. 57 Issue (8): 2644-2655   PDF    
基于照明补偿的单程波最小二乘偏移
周华敏, 陈生昌, 任浩然, 王汉闯, 陈国新    
浙江大学地球科学系, 杭州 310027
摘要:最小二乘偏移是一种基于反射地震数据与地下反射率间线性关系而建立起来的地震数据线性反演方法,相比常规偏移成像具有更好的保幅性能.本文提出了一种基于照明补偿的单程波最小二乘偏移方法,首先利用单程波方程的稳定Born近似广义屏波场传播算子构建反射地震数据与地下反射率间的线性算子,然后再应用线性最优化方法求解最小二乘偏移所对应的线性反问题.在迭代求解最优化问题的过程中,以地震波场的地下照明强度作为迭代反演的预条件算子加快迭代的收敛速度.单程波传播过程中考虑了速度分界面产生的透射效应,并用单极震源代替常规偏移中的偶极震源.把本文提出的方法应用于层状理论模型和Marmosi模型地震数据的数值试验中均取得了理想的结果.
关键词单程波     最小二乘偏移     照明强度     广义屏     振幅均衡性    
One-way wave equation least-squares migration based on illumination compensation
ZHOU Hua-Min, CHEN Sheng-Chang, REN Hao-Ran, WANG Han-Chuang, CHEN Guo-Xin    
Department of Earth Sciences, Zhejiang University, Hangzhou 310027, China
Abstract: Least-squares migration is a linear seismic inversion method which is established based on the linear relationship between seismic reflection data and underground reflection. Compared with conventional imaging technique, it has better performance in preserving amplitude. This paper presents a one-way wave equation least-squares migration based on illumination compensation. Firstly, we use the stable Born approximation generalized screen one-way wave-field propagator to construct the linear operator between reflected seismic data and underground reflection. Then we use the linear optimization method to solve the linear inverse problem of the least-squares migration. In the process of solving the optimization problems iteratively, we use the seismic wave illumination as a preconditioned operator of the iterative inversion, so as to accelerate the convergence speed of iteration. During the one-way wave propagation, we consider the transmission effect arising from velocity interface, and substitute monopole source for dipole source in conventional migration. Applying the method proposed in this paper to the layered theoretical model and the Marmosi model has achieved ideal results.
Key words: One-way wave     Least-squares migration     Intensity of illumination     Generalized Screen Propagator     Amplitude balance    

1 引言

油气储层附近反射体的精确成像对油气勘探有着至关重要的作用.但目前常规的偏移用正演算子的伴随代替其逆算子(Claerbout,1992),使得成像结果受观测系统、地下地质结构和地震数据的影响(Nemeth et al., 1999; Duquet et al., 2000)只能反映地下反射界面的位置,而不能获得准确的反射体振幅信息.为了使成像结果的保幅性更好,较为直接的办法是修改偏移算子以进一步接近正演算子的逆,这就是近年来发展较快的真振幅偏移成像技术(Etgen et al., 2009).该技术着重于研究真振幅的传播算子,通过补偿透射损失、几何扩散等来提高模拟算子的精度,因此该方法本质上没有突破常规偏移方法的局限性,不能达到真正保幅偏移的目的.Wu等(2004)提出了一种在角度域利用采集孔径效应进行偏移振幅校正的方法,陈生昌等(2007)提出了一种偏移阴影的方向照明补偿方法.波场照明补偿理论(Xie et al., 2006; 陈生昌和王汉闯,2010)对增强深层介质能量,提高成像振幅的均衡性有一定帮助,但这种局部角度域的偏移成像阴影补偿方法计算复杂,且单次偏移的照明补偿对于成像结果的改善有限.Gray(1997)比较了Delft方法、CWP方法和最小二乘偏移这三种振幅保真的偏移方法,并指出最小二乘偏移是保幅性能较高的一种方法.

最小二乘偏移(Least-Squares Migration,LSM)是在线性反演理论的指导下,应用线性最优化方法和迭代的方式求解最小二乘偏移所对应的线性反问题,以达到提高成像振幅保真度和成像分辨率的目的.Tarantola(1984)提出构建最小二乘的误差函数来求解反演问题.早期最小二乘偏移主要集中于发 展最小二乘Kirchhoff偏移.LeBras和Clayton(1988)Lambara等(1992)Dequet等(2000)Fomel等(2002)分别基于射线理论对最小二乘偏移进行了不同的探讨,Nemeth(1999)对最小二乘Kirchhoff偏移方法进行了系统的研究和总结,黄建平等(2013)基于最小二乘Kirchhoff偏移研究了碳酸盐岩裂缝 型储层的成像问题,刘玉金等(2013)讨论了局部倾角约束的最小二乘Kirchhoff 偏移.最小二乘Kirchhoff偏移能减少偏移假象,但Kirchhoff波场传播算子的计算精度较低,不能满足实际生产的需求.Kuehl和Sacchi(199920012003)提出了将最小二乘偏移应用到波动方程波场传播算子中,并将基于波动方程的最小二乘偏移用于AVA反演.Rickett(2003)通过傅里叶有限差分单程波最小二乘偏移来探讨地下的照明问题.沈雄君和刘能超(2012)详细讨论了裂步法最小二乘偏移的实现过程,杨其强和张叔伦(2008)实现了最小二乘有限差分法偏移.对于精度更高的广义屏单程波的最小二乘偏移仍未见有研究.近年来,基于全波方程的最小二乘偏移逐渐成为研究的热点(Dong et al., 2012).Dai等(2010)讨论了最小二乘逆时偏移中的若干实现问题.为了提高计算效率,Tang和Biondi(2009)提出了面向目标体的混合震源最小二乘偏移.虽然全波方程的最小二乘偏移精度较高,但计算效率低和对初始速度模型的依赖性强限制了其在最小二乘偏移中的应用;而波动方程的单程波传播算子计算效率较高,精度也能满足实际生产的需求.因此在当前的计算能力下,权衡计算效率与计算精度,采用单程波方程实现波动方程最小二乘偏移是一种比较实际的选择.

最小二乘偏移研究的难点之一在于求取Hessian算子及其逆(Chavent and Plessix, 1999; Lecomte,2008).目前多采用近似的方法求解(Duquet et al., 2000),但如何选取近似程度更高的Hessian算子近似值仍未有统一的意见(Plessix and Mulder, 2004; Valenciano et al., 2006; Ren et al., 2011).陈生昌等(2007)提出的照明补偿理论实际上就是考虑Hessian算子的作用,并以地下的照明分布近似代替Hessian算子,但该方法并未进行最小二乘偏移的迭代运算,因此对成像精度的提高有限.

在前人研究的基础上,本文基于地震数据的波动方程线性反演理论,提出了基于照明预条件的广 义屏最小二乘偏移方法(Preconditioned Generalized Screen Propagator Least-Square Migration,PGSP-LSM).以单程波方程的稳定Born近似广义屏波场传播算子作为最小二乘偏移方法中的波场传播算子,用地震波的地下双向照明强度作为最小二乘偏移中Hessian矩阵的近似,并把它作为最小二乘迭代求解过程中的预条件算子以改善收敛速度;再结合目前真振幅偏移的思想引入WKBJ近似对广义屏波场传播算子进行修改,考虑波场穿过速度分界面时波的透射效应,用单极震源代替常规偏移中的偶极震源.在理论层状模型和复杂Marmosi模型上的试验表明了本文提出的基于照明预条件的广义屏最小二乘偏移方法的有效性. 2 PGSP-LSM的方法原理 2.1 LSM的基本理论框架

基于地震波场的线性Born近似,反射地震数据的正问题可表示为(不考虑噪声):

其中d(xs,xg,ω)为该观测系统记录的地震数据,f(ω)为震源子波,G(x0,xs,ω)为从震源xs到反射点x0的格林函数,G(xr,x0,ω)为从反射点x0到检波点xr的格林函数,m(x0)为地层反射系数,Ω为成像区间.式(1)也可表述为:

其中,L =f(ω)∫ΩG(xr,x0,ω)G(x0,xs,ω)dV0是与采集形式、震源子波和速度密度模型有关的线性正传播算子,m 是地下反射率模型,d 为记录到的地震数据.

通过构建最小二乘意义下的误差泛函:

可以得到模型 m 的最小二乘解估计

式(4)所对应的解的估计也称为最小二乘偏移成像的结果. L T是 L 的伴随算子,对于我们要考虑的地震数据偏移成像问题,L和L T分别对应着地震波传播的模拟和偏移算子. L T L 就是Hessian算子,也 称为照明算子、去模糊化算子等(Aoki and Schuster, 2009; Schuster and Hu, 2000).由于 L T L 的规模较大,目前的计算条件下难以求得该算子及其逆的准确值,因此常采用迭代的方法得到式(4)的最小二乘解估计.目前比较常用的迭代方法是L and weber迭代法,其迭代格式为:

其中,k为迭代次数;mk为第k次迭代结果,mk+1为第k+1次迭代结果,αk为第k次迭代的步长.在迭代中,通常取m0=0.

L and weber迭代法是求解线性反问题的一般化公式,通过对步长α采取不同的计算策略,可以得到基于不同迭代格式的最小二乘偏移最优化解.若α取第k次迭代的Hessian的逆H-1k,则式(5)对应的是基于牛顿法的迭代格式:

式中,H = L T L 对应的就是Hessian算子,即反演误差泛函(3)式对模型 m 的二阶导数.

当取α为单位矩阵 I,则式(5)对应的就是基于梯度法的迭代格式:

如果对(7)式不进行迭代,得到的就是常规偏移结果:

由以上可知,常规偏移成像用正传播算子 L 的伴随算子 L T来代替逆算子(L T L)-1 L T,避免了算子求逆的过程,提高了求解过程的稳定性.但是受观测系统和上覆介质的影响,缺少(L T L)-1项的常规偏 移会导致成像分辨率低、振幅保幅性差.如果在式(6)中,用地下照明作为H-1k的近似来对常规偏移结果进 行修正,则得到照明补偿的偏移结果(陈生昌等,2007). 2.2 广义屏(GSP)传播算子的偏移与反偏移

偏移和反偏移是实现最小二乘偏移关键的两步,因此波场传播算子的精度对最小二乘偏移的成像精度有着至关重要的影响.目前主要是在双域中构造单程波传播算子,广义屏(GSP)波场传播算子(Wu and de Hoop,1996)相比裂步傅里叶(SSF)传播算子(Stoffa et al., 1990),可以实现波场在速度强横向变化的介质中比较精确的大角度传播,对介质的空间变化具有很好的适应性,而且计算的复杂性和计算量相比于傅里叶有限差分传播算子(FFD)(Ristow and Ruhl, 1994)更低.为此,本文选择采用稳定的Born近似广义屏波场传播算子来实现单程波最小二乘偏移(陈生昌,2002).

为了把zi深度层位上的波场p(x,zi;ω)传播到zi+Δz深度层上得到波场p(x,zi+Δz;ω)(Δz是传播间隔),可以利用已知波场p(x,zi;ω)来计算zi+Δz上的背景波场p0(x,zi+Δz;ω)和散射波场ps(x,zi+Δz;ω),其中背景波场p0(x,zi+Δz;ω)由下述相移方程给出:

式中,kz=为垂直波数,v0(zi)表示介质在深度间隔(zi,zi+Δz)内的背景参考速度,通常取平均速度,kx为水平波数,ω为频率.Fx和F<sup>-1<sub>kx分别代表x方向的傅里叶正变换和反变换.

利用局部Born近似并考虑波场穿过速度分界面时,界面对波场产生的透射效应,引入 WKBJ近似(Clayton et al., 1981),得到散射波场:

其中,ε(x,zi)=1/v2(x,zi)-1/v20(x,zi),k0(zi)=ω/v0(zi),a=k0(zi)/ki(zi).再利用数学近似式eiω v0(zi)ε/2(x,zi)Δz-1近似代替iω v0(zi)ε/2(x,zi)Δz,得到稳定的Born近似总波场计算公式:

p(x,zi+Δz;ω)=F-1kx eikz(zi)Δz a×Fx eiω v0(zi)2 ε(x,zi)Δz-1 p(x,zi;ω)+Fx p(x,zi;ω).(11)

这就是稳定的Born近似广义屏波场正传播算子 L .反传播算子 L T是正传播算子 L 的共轭转置.试验验证,稳定的Born近似广义屏波场传播算子在速度强横向变化介质的偏移成像中能取得很好的成像结果.

为了进一步提高最小二乘偏移的振幅保真性,这里对上述单程波传播算子进行修改.在张关泉保幅解耦思想指导下采用张宇(2006)刘定进(2007)等人提出的透射效应校正方法实现波场延拓过程中的振幅补偿功能,可以提高传播算子的保幅性.对波场在频率-波数域中做下述校正:

根据Wapenaar(1990)提出的单程波方程中震源表示,把常规偏移中的偶极震源变换成为单极震源,以进一步提高偏移的保幅性.

其中,S(ω)为震源子波的频谱.

根据反射地震数据的正演表达式(1)和地震数据的深度偏移过程,利用上文的稳定Born近似广义屏单程波传播算子进行地震数据的正演模拟,该过程也称为反偏移.首先把震源产生的入射波场向下传播,并记录每个深度网格点的入射波场值,直到最大深度层,再将最大深度处的入射波场与该层的反射率作用得到反射波场,然后向上正向传播反射波场,把每一层的入射波场与当前层的反射率作用得到的当前层的反射波场与正向传播的反射波场累加,直到到达地面检波点接收得到反射波场.对于广义屏单程波正演模拟时易出现的倏逝波、边界反射等问题,采取了相应的解决措施.采用广义屏算子的单程波正演模拟仅考虑地震波传播过程中的上行波和下行波中的一次反射波和不规则点的绕射波,算法简单快速,得到的地震记录信噪比高. 2.3 基于照明预条件的广义屏最小二乘偏移

为了提高迭代收敛速度,广义屏最小二乘叠前深度偏移可采用高斯-牛顿法的迭代格式,即先求取Hessian的一个不甚精确的近似值作为最小二乘偏移实现过程中的预条件,对迭代更新的成像结果做加权运算,从而达到提高成像精度、加快迭代收敛速度的目的,这就是预条件的最小二乘偏移.以地下的照明强度作为一个预条件算子来近似代替Hessian算子是一个很好的选择(陈生昌等,2007),即以炮检结合的照明结果作为最小二乘保幅偏移中Hessian矩阵的近似代替:

根据 L =f(ω)∫ΩG(xr,x0,ω)G(x0,xs,ω)dV0,可以得到:

为了提高计算效率,这里只考虑与地震数据主频对应的地震波的照明强度.在照明阴影区,利用照明分布对偏移结果进行照明补偿容易导致成像振幅的不稳定,产生成像噪声.为了提高后续求逆运算的稳定性,我们对(15)式增加一个数值极小的稳定因子ε,由此得到Hessian算子的近似计算表达式:

ε是一个数值极小的稳定因子,也称为阻尼因 子,ε的选取是否合适会直接影响成像结果,取值偏 大可能会改变成像振幅的比例关系,偏小又不能起到稳定作用.经过多次试验得到,当ε值取为地下照明强度最大值的0.001倍时,比较能符合数值计算的要求.由Hessian算子的主对角线元素可知,地下的照明强度实际上就是Hessian矩阵的主对角线元素,它综合了地震波场正传播和反传播的所有信息,能够反映地震观测系统对地下介质的照明情况以及对模型的分辨能力,因此也常被称为分辨率矩阵.本文的基于照明预条件的广义屏最小二乘偏移技术(PGSP-LSM)采用的就是介于常规的最速下降法与牛顿法之间的优化算法,一方面它可以避免直接求取Hessian算子的复杂性,另一方面通过照明预处理来加快LSM的收敛速度,提高成像振幅的均衡性,得到横向分辨率更高的反演成像结果.由此可得到PGSP-LSM的算法实现流程图:
图 1 基于照明预条件的广义屏最小二乘偏移(PGSP-LSM)的算法流程图 Fig. 1 The algorithm flow chart of the illumination Preconditioned Generalized Screen Propagation Least Squares Migration(PGSP-LSM)

图 1所示的流程图中可以看出,以最小二乘为基础的保幅偏移方法本质上等价于一种已知背景速度模型的线性反演方法.首先通过波场传播算子对炮检点的波场向下延拓,利用互相关成像条件计算地下成像点的成像结果,估求一个初始的地下反射率模型;再利用这个反射率模型重构地震数据,与实际地震记录进行匹配后,对数据残差进行成像;最后利用数据残差的成像结果对初始的反射率模型进行更新,对这一过程不断地迭代以后得到最优化的地下反射率估计.

地下的照明强度只与已知的背景速度模型有关而与迭代次数无关,即在迭代的过程中均采用同样的照明强度作为Hessian的近似值来对迭代的梯度方向进行预处理.单程波方程最小二乘偏移在进行数据残差的求取之前,应该做好单程波重构的地震数据与实际地震数据的匹配.由于常规的偏移成像中缺少Hessian算子项,而且重构的地震数据采用的广义屏单程波传播算子与实际地震数据也会有一定的差异.为了数据残差的求取,必须要将重构数据匹配到与实际数据同一个量级上.本文采取的措施是求取实际数据的均方根值与单程波正演模拟数据的均方根值的比值作为匹配因子来修正模拟数据.

综上所述,本文提出了采用基于照明预条件的广义屏最小二乘偏移(PGSP-LSM)方法来提高成像振幅的均衡性和保真性.考虑广义屏算子穿过速度分界面时的透射效应,并以单极震源代替偶极震源,用地下的双向照明结果近似最小二乘偏移成像中Hessian矩阵以提高反演迭代的收敛速度,从而达到增强深部目标体的波场能量强度、提高成像振幅均衡性的目的,理论上可以使成像结果具有更高的成像精度和成像分辨率. 3 数值试验

为了验证本文提出的基于照明预条件的广义屏最小二乘偏移方法的正确性和有效性,我们用一个简单的均匀层状模型和复杂的Marmousi模型进行数值试验. 3.1 简单层状模型试验

设计一个简单的横向均匀层状模型来验证本文方法的有效性.层状的速度模型如图 2a所示,模型大小为12 km×4 km,水平、垂直网格间距都是10 m.观测系统为中间放炮两边接收,共有300炮,起始的炮点位置在3 km,炮间距20 m,炮点的分布范围在3~9 km.每一炮的总道数为501道,道间距10 m,最小偏移距为0 m,记录长度3 s,时间采样率1 ms,记录数据是用声波方程的全波模拟得到的.图 2b为该层状模型的总照明强度图,计算照明强度采用的频率为15 Hz.

图 2(a)层状速度模型;(b)层状模型的总照明强度图 Fig. 2(a)Layered velocity model;(b)Total illumination of layered model

图 2b可以看出,地下的照明强度随深度的增加而减弱,浅层的照明最强,最底层能量最弱.照明强度在水平方向上的分布范围随着深度的增加而变大,而由于观测覆盖次数的减少,在观测系统两侧的照明强度没有中间覆盖次数较多的部分能量强,同一个深度层上越往两侧能量越弱. 图 3a为用声波方程进行全波模拟得到的地震记录,图 3b为采用本文的广义屏单程波传播算子模拟得到的正演记录.可以看出,本文的单程波传播算子具有较高的精度,除了少量能量较弱的边界反射外,已基本消除边界反射,在相位上能够很好地与基于双程波方程的模拟数据相吻合.为了进一步验证本文的单程波正演方法的有效性,对图 3的两种正演记录分别取出偏移距为2 km处的单道地震记录的波形振幅图进行对比(如图 4所示).图 4b的单程波正演的单道波形在相位上基本与图 4a中的声波模拟相一致,但是仍可以看出单程波模拟的记录波形在振幅值上与声波模拟并没有完全一致,并且存在一些小的噪声,这是单程波方法数值模拟所固有的缺陷,可以通过上文所述的数值匹配工作来尽量减少两者的振幅差异.由此说明采用本文的广义屏单程波传播算子进行最小二乘偏移是能够满足精度要求的.

图 3(a)单炮全波模拟地震记录;(b)本文所采用的广义屏传播算子的单程波模拟记录 Fig. 3(a)Full wave modeling seismic data;(b)Modeling seismic data of generalized screen propagation used in this paper

图 4(a)声波模拟地震记录;(b)单程波正演模拟地震记录在偏移距为2km处的波形振幅图 Fig. 4(a)Simulated seismic record of acoustic waves;(b)Waveforms of one-way wave simulated seismic record at 2 km offset

应用本文的基于照明补偿预条件的保幅最小二乘偏移(PGSP-LSM)方法于图 2a所示的层状模型中,得到常规的偏移结果和PGSP-LSM迭代5次的成像结果,截取成像图中水平坐标从2 km到10 km之间的成像结果显示在图 5中.对比图 5的两幅图可以看出,最小二乘偏移经过5次迭代的成像结果较常规偏移成像的振幅均衡性有了明显的提高,成像图中深层反射体的能量明显增强.在水平方向上,振幅的均衡性更强了,特别如图 5b中所圈出的部分,由常规偏移得到成像图中的振幅能量较弱,越往模型的两端,成像的振幅越小,甚至消失.对于PGSP-LSM的成像结果,在模型的两端,成像的振幅仍能保持在与中间位置相接近的振幅值,提高了非满覆盖区域的成像能量,使得成像能量分布有一定的横向均衡作用.此外,成像噪声有了明显的减少,模型的成像分辨率有了显著的提高.同时也可以看出在PGSP-LSM成像图 5b中沿偏移距方向,在非满覆盖次数的模型两端的浅层,由于照明强度较弱使得照明补偿后的成像振幅在小范围内出现不稳定现象.最小二乘偏移通过多次迭代能减少这种边界成像噪声,而基于计算效率的考虑,仅迭代5 次的成像结果中仍存在这种边界成像噪声,但这对成像效果影响不大.总体来说,本文的PGSP-LSM法能显著提高成像振幅的保真性、均衡性和成像分辨率.

为了仔细对比,我们抽取了图 5中各幅图相对应的6 km处的单道成像结果,如图 6所示.图 6a是由已知速度模型计算得到的理论反射率模型,图 6b是常规偏移成像结果,图 6c是最小二乘偏移迭代5次后的结果.由图 6的对比,可以更清楚地看到最小二乘偏移成像方法确实可以使深部偏移成像波场的能量得到更进一步的增强,在垂向上振幅的均衡性、保幅性得到了很大的提高.

图 5(a)常规偏移成像结果图;(b)迭代5次的PGSP-LSM成像结果 Fig. 5(a)Conventional migration imaging result;(b)PGSP-LSM imaging result after five iterations

图 6 图 5中单道的对比 (a)理论的反射率;(b)常规偏移结果;(c)PGSP-LSM迭代5次的结果. Fig. 6 Comparison of single channel from Fig. 5 (a)Theoretical reflectivity;(b)Conventional migration result;(c)PGSP-LSM result after 5 iterations.

为了对比各种方法在振幅保真性方面的效果,对图 5中的常规偏移成像结果和PGSP-LSM成像结果沿层拾取各个反射面上的最大振幅值,截取炮点覆盖的范围3~9 km的结果展示在图 7图 7a为三个反射面的理论反射率值,图 7b为采用常规偏移方法拾取的三个反射面的反射率值,图 7c为采用本文的PGSP-LSM方法拾取的三个反射面的反射率值.通过对比,可以很明显地看出常规偏移的保幅性较差,仅在观测覆盖次数较多的区域得到相对保幅的振幅值,在水平反射层的两侧则表现为明显的振幅不足,而且越往深层,反射体成像振幅值的估计明显低于实际值.而本文所采用的PGSP-LSM的振幅 保真性较常规偏移有了很大的提高,在垂直方向上表现为深层的反射面的振幅值有明显的增强,垂向振幅的均衡性有明显提高,沿反射层方向的模型两侧,成像振幅不足的问题得到显著改善,成像振幅值沿炮检距方向的变化趋势大体能反映地下反射界面反射率的变化情况.选取某一个反射界面(以第二层反射面为例)的上述几种成像方法得到的反射率的对比(如图 7d所示),也可以看出本文的PGSP-LSM方法在水平方向的振幅延续度较常规偏移更好,在一个更大的水平偏移距范围内振幅都能保持在一个相似的水平上,这与地下实际反射层的反射率变换情况相吻合,振幅在横向上的保真度也得到了很好的增强.但是也可以从图 7c中看出,由于本项研究采用近似程度不高的Hessian近似值以及存在边界效应,使得基于照明预条件的最小二乘偏移成像结果两侧的照明阴影区出现了成像振幅的不稳定现象,具体表现为浅层是振幅欠补偿,深层是振幅过补偿.这个问题需进一步通过采用精度更高的Hessian近似值、去除边界效应等措施来解决.

图 7 反射率估计试验对比 (a)理论反射率;(b)常规偏移;(c)本文的PGSP-LSM方法;(d)选取第二层对比几种方法的反射率. Fig. 7 Comparison of different methods in estimating reflectance (a)Theoretical reflection;(b)Conventional migration;(c)PGSP-LSM method used in this paper;(d)Second reflection estimated by three methods.

通过该简单的层状模型可以证明,最小二乘偏移确实对增强深部目标体的偏移成像振幅能量有很明显的效果,可用于解决深部目标体的偏移能量弱的问题,有助于提高成像振幅的保真度,从而实现真振幅偏移. 3.2 复杂模型试验

为了验证本文的基于照明预条件的最小二乘偏移方法的正确性和有效性,我们采用国际通用标准模型——Marmousi模型数据来进行数值试验.Marmousi模型的网格大小为343×750,水平网格间距为25 m,垂直网格间距为4 m,Marmousi模型数据共有240炮,每一炮有96个接收道,炮间距为25 m,道间距为25 m,首炮地面位置为2575 m,最小炮检距为200 m,最大炮检距为2575 m;记录长度为3 s,时间采样率为4 ms,主频为25 Hz.

应用2.3节提出的照明计算方法计算波场的照明强度,图 8b为Marmousi模型的地下照明分布.该照明图中表现出了较为复杂的照明不均匀性,但还是可以分辨出地下主要的三个断层构造.越往深层,照明能量就越弱.

图 8 Marmousi模型的(a)速度模型及(b)地下照明强度分布 Fig. 8 Velocity model(a) and illumination(b)of Marmousi model

应用本文提出的基于照明预条件的广义屏单程波最小二乘偏移方法对Marmousi模型进行成像.图 9a为Marmousi模型常规偏移的结果,图 9b为迭代10次的PGSP-LSM成像结果.本文的基于照明预条件的最小二乘偏移法只需要经过10次迭代得到的成像振幅的均衡性就能较常规偏移成像结果有明显的增强,浅层陡倾断层成像清晰,由于上覆浅层三大断层和中层背斜所引起的深层能量损失的现象有明显的改善,突显了深层局部构造,整体表现为能量更为均衡,有效改善成像结果质量,成像的分辨率也有明显的提高.图 9b的上部、中下部和右边的成像分辨率都有明显的提高.但也要注意在该图中的右下角出现了噪声放大的现象.为了仔细地对比,特抽取图 9的两幅图的第230道与该道的理论反射率进行对比,如图 10所示.

图 9 Marmousi模型的(a)常规偏移成像结果及(b)PGSP-LSM迭代10次成像结果 Fig. 9 Conventional migration result(a) and PGSP-LSM result after 10 iterations(b)of Marmousi

图 10 Marmousi模型(a)理论反射率,(b)常规偏移,(c)PGSP-LSM迭代10次成像的单道振幅对比 Fig. 10 Comparison of single channel’s amplitude for Marmousi model:(a)theoretical reflection, (b)conventional migration, and (c)PGSP-LSM result after 10 iterations

图 10为选取Marmousi模型的垂直反射率、常规偏移结果和本文所提的PGSP-LSM成像结果的单道记录的归一化振幅的对比图.可以更清楚地看到,最小二乘偏移在10次迭代后的深层反射体的能量有了明显提高,成像振幅均衡性、成像分辨率较常规偏移(图 10b)显著增强.而且与垂直反射率(图 10a)相比,最小二乘偏移得到的反射率是与入射角无关的平均反射率,其单道成像结果综合了地震波在各个方向上的传播信息,因此相比垂直反射率能更真实客观地反映地下介质的真实变化信息.图 10c中1.8 km以下的地下介质深层部分,常规偏移 不能很好地成像,而采用本文的PGSP-LSM方法则能显著提高成像的振幅,且与理论的反射率相对应.由浅到深,整体的振幅均衡性有明显的提高,在部分区域成像分辨率也提高了.但是从图 9也可以看出,最小二乘偏移在经过多次迭代后,增加了一些成像噪声,这是采用基于照明预条件的最小二乘偏移所不能 避免会产生的问题,这也是需要我们后续努力解决的.

和常规偏移相比较,最小二乘偏移会带来一定的计算量.一般来说,最小二乘偏移完成1次迭代运算的计算时间约为常规偏移的2倍.但是由以上两个例子可以看出,本文所提的基于照明预条件的最小二乘偏移方法只需要经过少量的迭代即能较常规偏移取得更为满意的成像结果.随着并行技术的发展、高性能计算机群的使用,最小二乘偏移所带来的计算效率问题可以得到极大的缓解.当前单程波方程偏移成像技术在实际生产单位中有着广泛的应用,而基于单程波方程的最小二乘偏移技术只需在常规偏移技术的基础上经过多次迭代即能获得更好的成像结果.尤其在常规偏移难以精确成像的复杂构造区,最小二乘偏移能够得到更为均衡的、反映地下界面信息的成像结果,这有利于将该技术推广应用于实际地震资料处理解释中.虽然本文的工作和研究结果都是基于二维形式,但是它们都能非常简便地推广到三维情况. 4 结论

本文提出的基于照明预条件的广义屏最小二乘偏移是在线性反演理论的指导下,采用稳定的Born近似广义屏算子作为单程波最小二乘偏移的传播算子,以炮检结合的照明结果作为最小二乘保幅偏移中Hessian算子的近似,考虑波场传播穿过分界面的透射效应,并用单极震源代替常规偏移中的偶极震源,通过迭代运算求解最小二乘意义下的最优化问题.在理论模型和复杂的国际标准模型上的试验证明了本文所提的最小二乘偏移方法在经过少量的迭代后,确实能较常规偏移成像提高成像振幅的保真性、均衡性,提高成像分辨率,可用于解决深层目标体成像能量弱的问题,有助于增强弱信号和提高复杂构造的成像质量.

尽管在理论模型上的试验取得了较满意的结果,但是本文方法存在的噪声问题是不可忽视的.随着迭代次数的增加,通过照明强度对成像结果进行预处理会带来一些累加噪声,这是该方法自身所不可避免的,如何能减少甚至消除这些噪声是我们下一步应该进行研究的内容.此外,如何提高最小二乘偏移的计算效率,并开展成像精度更高的基于全波方程的最小二乘偏移也是我们下一步应着重研究的方向.总的来说,最小二乘偏移作为传统的偏移成像与基于非线性反演理论的全波形反演之间的桥梁和纽带,对其进行深入研究有助于我们更深刻认识地震波传播与成像,最终实现对地下进行精确的成像.

参考文献
[1] Aoki N, Schuster G T. 2009. Fast least-squares migration with a deblurring filter.   Geophysics, 74(6): WCA83-WCA93.
[2] Chavent G, Plessix R E. 1999. An optimal true-amplitude least-squares prestack depth-migration operator.   Geophysics, 64(2): 508-515.
[3] Chen S C. 2002. Research on methods for wave equation pre-stack depth migration and velocity model construction (in Chinese). Shanghai: Tongji University.
[4] Chen S C, Ma Z T, Wu R S. 2007. Mumination compensation for wave equation migration shadow. Chinese J. Geophys. (in Chinese), 50(3): 844-850.
[5] Chen S C, Ma Z T, Wu R S. 2007. Two-way subsurface directional illumination analysis by wave equation.   Journal of Tongji University: Natural Science (in Chinese), 35(5): 681-684.
[6] Chen S C, Wang H C. 2010. Migration compensation with plane wave illumination. Chinese J. Geophys.   (in Chinese), 53(7): 1710-1715.
[7] Claerbout J F. 1992. Earth Soundings Analysis: Processing Versus Inversion. Blackwell Science.
[8] Clayton R W, Stolt R H. 1981. A Born-WKBJ inversion method for acoustic reflection data.   Geophysics, 46(11): 1559-1567.
[9] Dai W, Boonyasiriwat C, Schuster G T. 2010. 3D multi-source least-squares reverse time migration.  //Proceedings of the 2010 SEG Annual Meeting, 29: 3120-3124.
[10] Dong S, Cai J, Guo M, et al. 2012. Least-squares reverse time migration: towards true amplitude imaging and improving the resolution.  //SEG Technical Program Expanded Abstracts 2012, Society of Exploration Geophysicists, 1-5.
[11] Duquet B, Marfurt K J, Dellinger J A. 2000. Kirchhoff modeling, inversion for reflectivity, and subsurface illumination.   Geophysics, 65(4): 1195-1209.
[12] Etgen J, Gray S H, Zhang Y. 2009. An overview of depth imaging in exploration geophysics.   Geophysics, 74(6): WCA5-WCA17.
[13] Fomel S, Berryman J G, Clapp R G, et al. 2002. Iterative resolution estimation in least-squares Kirchhoff migration.   Geophysical Prospecting, 50(6): 577-588.
[14] Gray S H. 1997. True-amplitude seismic migration: A comparison of three approaches.   Geophysics, 62(3): 929-936.
[15] Huang J P, Li Z C, Kong X, et al. 2013. A study of least-squares migration imaging method for fractured-type carbonate reservoir. Chinese J. Geophys. (in Chinese), 56(5): 1716-1725.
[16] Kuehl H, Sacchi M D. 1999. Least-squares split-step migration using the Hartley transform. Univ.   of Alberta, SEG Expanded Abstracts, 1548-1551.
[17] Kuehl H, Sacchi M F. 2001. Split-step WKBJ least-squares migration/inversion of incomplete data. Proceedings of the 5th SEG International Symposium-Imaging Technology.
[18] Kuhl H, Sacchi M D. 2003. Least-squares wave-equation migration for AVP/AVA inversion.   Geophysics, 68(1): 262-273.
[19] Lambare G, Virieux J, Madariaga R, et al. 1992. Iterative asymptotic inversion in the acoustic approximation.   Geophysics, 57(9): 1138-1154.
[20] LeBras R, Clayton R W. 1988. An iterative inversion of back-scattered acoustic waves.   Geophysics, 53(4): 501-508.
[21] Lecomte I. 2008. Resolution and illumination analyses in PSDM: A ray-based approach.   The Leading Edge, 27(5): 650-663.
[22] Liu D J. 2007. Preserved-amplitude seismic migration imaging method based on wave equation (in Chinese). Beijing: China University of Petroleum.
[23] Liu Y J, Li Z C, Wu D, et al. 2013. The reserch on local slope constrained least-squares migration. Chinese J. Geophys. (in Chinese), 56(3): 1003-1011.
[24] Nemeth T, Wu C, Schuster G T. 1999. Least-squares migration of incomplete reflection data.   Geophysics, 64(1): 208-221.
[25] Plessix R E, Mulder W A. 2004. Frequency-domain finite-difference amplitude-preserving migration.   Geophysical Journal International, 157(3): 975-987.
[26] Ren H, Wu R S, Wang H Z. 2011. Wave equation least square imaging using the local angular Hessian for amplitude correction.   Geophysical Prospecting, 59(4): 651-661.
[27] Rickett J E. 2003. Illumination-based normalization for wave-equation depth migration.   Geophysics, 68(4): 1371-1379.
[28] Ristow D, Rühl T. 1994. Fourier finite-difference migration.   Geophysics, 59(12): 1882-1893.
[29] Schuster G T, Hu J X. 2000. Green's function for migration: Continuous recording geometry.   Geophysics, 65(1): 167-175.
[30] Shen X J, Liu N C. 2012. Split-step least-squares migration.   Progress in Geophysics (in Chinese), 27(2): 761-770.
[31] Stoffa P L, Fokkema J T, de Luna Freire R M, et al. 1990. Split-step Fourier migration.   Geophysics, 55(4): 410-421.
[32] Tang Y, Biondi B. 2009. Least-squares migration/inversion of blended data. Proceedings of the SEG Technical Program Expanded Abstracts, 2859-2863.
[33] Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation.   Geophysics, 49(8): 1259-1266.
[34] Valenciano A A, Biondi B, Guitton A. 2006. Target-oriented wave-equation inversion.   Geophysics, 71(4): A35-A38.
[35] Wapenaar C P A. 1990. Representation of seismic sources in the one-way wave equations.   Geophysics, 55(6): 786-790.
[36] Wu R S, Chen S C, Luo M Q. 2004. Migration Amplitude Corrections in Angle-Domain Using Beamlet Decomposition.   Proceedings of the 66th EAGE Conference & Exhibition.
[37] Wu R S, de Hoop M V. 1996. Accuracy analysis and numerical tests of screen propagators for wave extrapolation.   Proceedings of the SPIE's 1996 International Symposium on Optical Science, Engineering, and Instrumentation, International Society for Optics and Photonics, 196-209.
[38] Xie X B, Jin S W, Wu R S. 2006. Wave-equation-based seismic illumination analysis.   Geophysics, 71(5): S169-S177.
[39] Yang Q Q, Zhang S L. 2008. Least-squares Fourier finite-difference migration. Progress in Geophysics (in Chinese), 23(2): 433-437.
[40] Zhang Y. 2006. The theory of true amplitude one-way wave equation migration. Chinese J. Geophys.   (in Chinese), 49(5): 1410-1430.
[41] 陈生昌. 2002. 波动方程叠前深度偏移成像与速度建模方法研究[博士论文].   上海: 同济大学.
[42] 陈生昌, 马在田, Wu Ru Shan. 2007. 波动方程偏移成像阴影的照明补偿.   地球物理学报, 50(3): 844-850.
[43] 陈生昌, 马在田, 吴如山. 2007. 波动方程双程地下方向照明分析.   同济大学学报: 自然科学版, 35(5): 681-684.
[44] 陈生昌, 王汉闯. 2010. 基于平面波照明的偏移成像补偿.   地球物理学报, 53(7): 1710-1715.
[45] 黄建平, 李振春, 孔雪等. 2013. 碳酸盐岩裂缝型储层最小二乘偏移成像方法研究.   地球物理学报, 56(5): 1716-1725.
[46] 刘定进. 2007. 波动方程保幅地震偏移成像方法研究[博士论文].   北京: 中国石油大学.
[47] 刘玉金, 李振春, 吴丹等. 2013. 局部倾角约束最小二乘偏移方法研究.   地球物理学报, 56(3): 1003-1011.
[48] 沈雄君, 刘能超. 2012. 裂步法最小二乘偏移.   地球物理学进展, 27(2): 761-770.
[49] 杨其强, 张叔伦. 2008. 最小二乘傅里叶有限差分偏移. 地球物理学进展, 23(2): 433-437.
[50] 张宇. 2006. 振幅保真的单程波方程偏移理论.   地球物理学报, 49(5): 1410-1430.