2016, Vol. 59
2. 中海油研究总院, 北京 100028
2. CNOOC Research Institute, Beijing 100028, China
VSP 技术在井旁构造精细成像、井旁断层识别、井旁地层岩性描述和地震波衰减等方面都取得了许多实际的应用效果.VSP成像是VSP 数据处理中的一个重要环节，VSP成像方法主要包括常规的VSP-CDP转换、基于射线理论的Kirchhoff积分法偏移和基于波动方程的偏移三种，前两种方法都有其局限性，VSP-CDP转换仅对水平层状介质或横向变速缓慢的介质才能取得比较好的处理结果，Kirchhoff 积分方法因它固有的高频近似假设和射线理论在复杂介质中的缺陷(如焦散和多路径等) 使其在复杂构造区难以取得理想的结果.而波动方程偏移成像方法相对于射线追踪技术，能更准确地描述波在复杂介质中的传播以及波场能量的变化情况，这使得基于波动方程的偏移方法更受青睐.波动方程偏移方法又可以进一步分为基于单程波方程的偏移方法和基于双程波方程的偏移方法.基于双程波的逆时偏移不受地下界面倾角的限制，并且在保幅方面具有优势.孙文博和孙赞东(2010)开展了VSP资料伪谱法逆时偏移研究，李文杰等(2012)研究了二维各向同性介质中两分量的非零偏VSP叠前逆时深度偏移方法，蔡晓慧等(2015)实现了基于自适应优化有限差分方法的全波VSP逆时偏移，但该方法主要存在的应用瓶颈是：计算量大、存储量大和噪声问题.在波场的正向和逆时外推过程中需要存储各个时刻的波场，非常耗时.同时将针对单程波方程提出的互相关成像条件用于逆时偏移，会产生低频噪音干扰.因此三维逆时偏移的计算量和对速度模型的高要求使得它在目前工业界的应用效果并不是十分理想.
Claerbout(1976)首次提出了时空域单程波方程，然而该方程具有明显的倾角限制，无法对高陡构造准确成像，一些学者(马在田，1983) 采用不同的高阶近似方法来提高成像倾角，但由于方程中往往包含高阶偏导数项，而高阶差分计算量明显增加，使得数值计算难以实现.Guddati(2006)提出了任意广角波动方程(AWWE)，它是一种在空间域表述的高精度单程波方程.AWWE的优点是精度高，要提高AWWE的成像倾角，只需增加参考速度的个数即可，AWWE的形式简单，只包含二阶偏导数项，易于数值实现，因此AWWE被广泛应用于偏移成像，何兵寿等(2008a)用AWWE对地面地震资料做逆时偏移成像研究，同时将其与基于双程声波方程的逆时偏移方法做比较，指出基于AWWE的偏移方法不仅可以对陡倾角的地层准确成像，而且还可以很好地避免双程波方程逆时偏移的低频成像噪音；孙歧峰和杜启振(2011)推导了频率-空间域的AWWE，并利用有限差分高阶分裂法求解方程，实现地震偏移成像，指出该方法能够适应强横向变速地区和陡倾角地层的成像.在吸收边界条件方面，Heidari和Guddati(2006)导出广角方程正向传播的吸收边界条件；何兵寿等(2008b)将加权校正及插值预测吸收边界条件引入任意广角方程的逆时偏移领域，这种边界条件的本质是对Clayton 吸收边界条件的改进，它难以有效地吸收地震波逆时延拓过程中投射到边界的外行波；何兵寿等(2010)给出了适用于二维任意广角单程声波方程逆时延拓的吸收边界条件，给出的算法只需在少量增加计算量的前提下有效地压制了截断边界的伪反射，消除了边界处的偏移噪声；Chen等(2013)将完全匹配层(PerfectlyMatched Layer，PML) 用于AWWE，有效地解决了AWWE数值计算中的吸收边界问题.尽管AWWE具有高精度的特点，但是其计算过程包含重复的矩阵求逆和矩阵相乘，使得其计算代价比其他单程波方程高出许多，特别是增加参考速度个数时，其计算量显著增加.周辉等(2014)改进了现有的有限差分计算方案，大幅度提高了计算效率，同时采用了一种有效的方法在f-k域比较准确地区分出倏逝波和非倏逝波，从而获得一个合理的视速度作为f-k滤波器的门槛来压制倏逝波干扰，提高成像质量.本文将AWWE推广应用至VSP领域，实现了VSP资料的叠前深度偏移，采用非线性算法提高方程的成像精度，并应用了PML吸收边界条件，很好地压制了边界反射成像噪音.
2 基本原理 2.1 高角度单程波VSP叠前深度偏移算法三维标量AWWE下行波方程的一般形式为(Guddati and Heidari，2005)
|
(1) |
其中 c 为背景速度， dT=(1，0，…，0)n×1，u=(u，u1，u2，…，un－1)T，u1，u2…，un－1 为辅助变量，T表示转置，x、y、z为空间坐标，t表示时间变量，
|
(2) |
|
(3) |
式中， c1，c2，…，cn 为根据实际介质速度选取的参考速度.
2.2 高角度单程波时域PML匹配方程在进行高角度单程波方程有限差分偏移时会遇到边界反射问题，因此针对下行波式(1) 推导其PML匹配方程，考虑沿x方向和y方向衰减，将其分裂成三项，即
|
(4) |
将式(4) 前两式变换到频率域，用
|
(5) |
其中 u1、u2、u分别为u1、u2、u 的Fourier变换，ω为角频率.式(5) 第一式两端同时乘以 (iω+dx(x))2， 第二式两端同时乘以 (iω+dy(y))2， 消去公共项 (iω)2， 得到
|
(6) |
返回时间域得到
|
(7) |
将式(7) 与式(4) 的第三式写在一起，即组成了三维任意广角单程波下行波的时域PML匹配方程，即
|
(8) |
由于是单程波方程不需要考虑Z方向的边界反射问题，在计算时只需考虑前后左右四个边界和四周四个棱边，X方向两个边界只需要考虑X方向衰减，Y方向两个边界只需要考虑Y方向衰减，而四个棱边需要同时考虑X和Y方向的衰减.
2.3 非线性反演算法提高方程成像精度在推导高角度单程波方程时，不可避免地要对平方根算子做近似，对平方根算子的近似程度直接影响偏移成像的精度.从平方根算子表达的频散关系出发，对参考速度的选择进行优化.优化的准则是在最小平方约束下使高角度方程对平方根算子的近似与准确的单程波平方根算子达到最佳的匹配.下面选定两个参考速度来详细说明高角度声波方程参考速度的优化过程.
下行单程波平方根算子的准确表达式可以写为
|
(9) |
令 
|
(10) |
其中θ表示波的传播方向与界面法线之间的夹角.
二维AWWE下行波方程：
|
(11) |
为了评价AWWE的精度，对上式中的x、t做傅里叶变换，就可以得到频率波数域中的高角度方程
|
(12) |
进而可以整理出高角度方程对平方根算子的近似函数 Fn(ξ，c1，…，cn)， 当选取两个参考速度，并对速度进行归一化后，得到 F2(ξ，c1，c2) 为(详细推导过程参见附录A)：
|
(13) |
(8) 与式(10) 的相对误差为
|
(14) |
在最小平方约束条件下，使得它们的相对误差 Er 最小.相对误差越小，平方根算子的近似程度越高，越有利于大角度地层准确成像.因此研究方程的精度对实际地震资料处理有一定的价值.可以通过非线性反演算法(牛顿法) 求得 Er 最小时的参考速度(参见附录B).
在两个参考速度的情况下，当固定第一个参考速度 c1=1 时，得到的优化参考速度是 c1=1，c2=5.24； 当不固定参考速度 c1 的数值时，得到的优化参考速度是 c1=1.23，c2=6.46. 在固定 c1=1 时的三个参考速度的情况下得到的优化参考速度是 c1=1，c2=2.82，c3=21.15. 图 1是四种参考速度下AWWE频散曲线，图中OWWE代表的是下行单程波算子准确表达式的频散曲线.图 1a是参考速度 c1=1，c2=4 的AWWE频散曲线，AWWE的精度可以达到75°左右.图 1b是参考速度为 c1=1，c2=5.24 的AWWE频散曲线，在80°以下它与准确的频散曲线基本上是符合的，与图 1a相比，图 1b的最大准确成像角度由75°提高到80°.图 1c是参考速度为 c1=1.23，c2=6.46 的AWWE频散曲线，与图 1a和图 1b相比，其精度可以达到82°左右，图 1d是参考速度为 c1=1，c2=2.82，c3=21.15 的AWWE频散曲线，最大准确成像角度可接近90°，同时可以看到AWWE的频散曲线与准确的频散曲线整体上都可以精确地吻合，成像精度大大得到了提高.
|
图 1 不同参考速度对应的AWWE频散曲线及其误差 (a)c1=1，c2=4;(b)c1=1，c2=5.24;(c)c1=1.23，c2=6.46;(d)c1=1，c2=2.82，c3=21.15，图中OWWE代表的是下行单程波算子准确表达式的频散曲线. Fig. 1 The dispersion relations and errors with different reference velocities |
从上面的分析中可以看出：优化的三个参考速度AWWE频散曲线的精度最高；优化的两个参考速度AWWE频散曲线的最大准确成像角度都得到了提高，但是会损失小角度处的精度；与Guddati给出的 c1=1，c2=4 参考速度相比，优化方法选取的参考速度可使大角度处AWWE的精度得到提高.采用优化参考速度参数的AWWE做叠前深度偏移处理，而非采用增加参考速度个数的方法，就可以在避免增加运算量的情况下进一步解决大角度成像的问题.
3 数值模拟实验 3.1 层状介质模型用一个四层介质模型(如图 2a所示) 的叠前深度偏移实验分析高角度单程波VSP数据的偏移效果.模型大小为1800 m×3000 m，界面分别在900 m、1800 m和2700 m深度处，四个水平层速度分别为2000 m·s-1 、2200 m·s-1 、2400 m·s-1和2600 m·s-1.空间网格大小为3 m，时间采样间隔为0.5 ms，合成地震记录长度为2 s.炮点深度为3 m，水平方向600个网格点上均匀放置600炮激发.井水平坐标903 m，从深度为300 m到深度1200 m，以15 m为间隔均匀布置60个检波点.
|
图 2 简单模型偏移结果 (a) 水平层状模型；(b) 偏移结果. Fig. 2 The migration result of simple model (a)Horizontal layered model；(b)Its migration result. |
我们利用二维声波方程正演得到原始的共炮集数据之后，由于直达波(下行波) 是先于反射波(上行波) 到达检波器的，所以我们可以在共炮道集上对每个道沿着时间增大的方向寻值，根据直达波信号的大小，当寻到第一个与直达波信号量级相当的数值时，沿着时间增大方向将一个子波长度内的值变为零，这样就可以去除直达波.去除了直达波之后抽取共检波点道集，再利用共检波点道集数据通过高角度单程波叠前深度偏移得到60个共检波点道集偏移结果，通过叠加得到最后的结果如图 2b所示，可以看出三个水平界面得到准确成像，偏移噪音很小，但是成像范围并不是很大，这是由VSP本身的观测范围所限制的.
3.2 复杂模型我们用复杂模型(如图 3a所示) 的叠前深度偏移实验分析高角度单程波VSP数据的偏移效果.模型大小为3000 m×3000 m，水平方向网格大小为5 m，深度方向网格大小为5 m，时间采样间隔为0.5 ms，合成地震记录长度为2 s.从上到下各地层中声波速度分别为4000 m·s-1、4100 m·s-1、4400 m·s-1、4100 m·s-1、4000 m·s-1、4200 m·s-1和4300 m·s-1.炮点深度为5 m，水平方向600个网格点上均匀放置600炮激发.井水平坐标2000 m，从深度为250 m到深度850 m，以20 m为间隔均匀布置30个检波点.
|
图 3 复杂模型偏移结果 (a) 复杂模型；(b) 偏移结果. Fig. 3 The migration result of complex model (a) Complex model；(b) Its migration result. |
图 3b是多炮偏移叠加成像的结果，可以看出，对于离井较近的高陡地层得到了很好地成像，偏移的成像噪音也得到了很好的压制.
4 实际资料测试选取Q油田Walk-away VSP测线进行了叠前深度偏移成像，本次Walk-away VSP 观测系统震源布设范围为-1900～1900 m，炮间距25 m，检波器布设范围为600～800 m，采用间距为10 m的20级检波器接收，每个“排列”上提5 m重复测量，最终达到5 m的接收间距，采集时一共布置了八条测线，这里对其中一条测线进行处理.图 4为偏移速度模型，由于只是对其中一条测线的处理，所以这里给出的只是截取了的一部分速度模型，并将坐标值均转化为正值，以便于描述.抽取共成像点道集进行偏移成像处理，同时为了更好地进行同相叠加，抽取了40个共检波点道集不同偏移距的共成像点道集，利用互相关方法对得到的共成像点道集进行整道校正，从而提高成像质量.图 5为水平偏移距400 m校正前后的共成像点道集，对比可以看出，校正后的共成像点道集同相轴连续性得到增强，图 6为共成像点校正后的偏移结果，可以看出同相轴连续性良好，分辨率也很高，这主要是得益于VSP数据具有较高的主频和信噪比.
|
图 4 偏移速度模型 Fig. 4 Velocity model for migration |
|
图 5 共成像点道集校正(偏移距400 m处) (a) 校正前；(b) 校正后. Fig. 5 The correction of CIP gathers(offset=400 m) (a) Before；(b) After. |
|
图 6 共成像点校正后的VSP偏移结果 Fig. 6 VSP migration result after CIP correction |
为了进一步分析成像结果，将地面地震偏移成像结果与此次VSP偏移结果进行对比，如图 7，可以看出VSP偏移结果与地面地震剖面产状吻合较好，并具有更高的分辨率和保幅性，主频达到100 Hz，频带宽度达到8~180 Hz，为后续波组特征的识别和断层的精细刻画提供了高品质资料基础.
|
图 7 VSP与地面地震偏移结果的对比 (a) VSP剖面与地面地震剖面叠合图；(b) 地面地震剖面频谱；(c) VSP剖面频谱. Fig. 7 The comparison of VSP migration result and surface seismic section (a) The overlapping graph of VSP migration result and surface seismic section；(b) Spectrum of surface seismic section；(c) Spectrum of VSP migration result. |
与其他单程波波动方程相比，AWWE是一种高精度的单程波方程，可以实现高陡构造的准确成像.本文将这种时空域高角度单程波方程偏移推广应用到VSP数据中，从三维标量任意广角波动方程出发，推导了PML公式，很好地压制了边界反射成像噪音，同时应用牛顿法在最小平方意义下求解方程组来选取参考速度，从而提高了地层成像倾角，优化后成像倾角接近90°.模型试算结果表明该方法在陡倾角情况下能取得较好的成像效果，实际资料的试处理也表明VSP时空域高角度单程波方程偏移成像结果相比地面地震偏移剖面具有更高的分辨率.
附录A 提高精度的理论分析二维时空域下行波AWWE为
|
(A1) |
AWWE作为一种单程波方程(OWWE)，不可避免地要对平方根算子做近似，其近似程度直接影响偏移成像的精度和所能准确成像的界面倾角范围.用水平波数kx，角频率ω表示的准确的下行单程波方程表达形式如下：
|
(A2) |
令 ckx/ω=sinθ，θ 表示波的传播方向与界面法线之间的夹角， θ 的范围是0到90°，所以(A2) 可以简化为
|
(A3) |
为了提高AWWE的精度，可将AWWE二维下行波方程化成上式 
所以我们对方程(A1) 中的x，t做傅里叶变换，就可以得到
|
(A4) |
其中 u是u关于x，t 的二维傅里叶变换，
|
(A5) |
移项可得
|
(A6) |
由于参考速度通常设置为两个，多个参考速度的情况参照此法也可同理得到优化参考速度.因此为了方便讨论，我们选择两个参考速度的情况，代入
|
|
(A7) |
AWWE的 Λ1与Λ2 的形式决定AWWE的偏移精度只是受参考速度与背景速度的比值影响，因此我们可以假定背景速度为常数(为方便讨论，本文假定背景速度为1).由上面矩阵方程可以得到两个方程，再通过消去u1 可以得到
|
(A8) |
相对误差 Er 最小时，平方根算子的近似函数 Fn(ξ，c1，…，cn) 精度最高. Er 最小时的参考速度即为我们要求取的优化参考速度.
|
(B1) |
定义目标函数
|
(B2) |
其中 θ1和θ2 是所选择的角度范围，通常选取θ1=0°，θ2=90°.目标函数Ir表示的是在一定角度范围内对误差 Er 的最小平方约束.目标函数越小，则F2与R越接近，此时AWWE的成像精度越高.
通过求取目标函数 Ir 的最小值，得到参考速度的值.同样以两个参考速度的情况为例，在最小平方意义下由方程(B2) 可得
|
(B3) |
方程组(B3) 是一个非线性方程组，可以用非线性反演算法(如梯度法、牛顿法、遗传算法、模拟退火算法等) 求解.本文选用牛顿法求解方程组(B3).首先给定c1和c2的初值c10和c20，然后采用方程组(B4) 更新参考速度.
|
(B4) |
其中 c1(k)，c2(k) 表示第 k 次的参考速度迭代值， Δc1(k)，Δc2(k) 表示第 k 次的参考速度更正值， c11(k+1)，c2(k+1)是第k+1 次的参考速度迭代值，参考速度的更正值 Δc1(k)，Δc2(k) 由方程(B5) 求出：
|
(B5) |
其中 Δc(k)=[Δc1(k)，Δc2(k)] 是参考速度更正向量， g(c(k))=g(c1)，g(c2)(k) 是目标函数的梯度向量，
|
(B6) |
|
(B7) |
|
(B8) |
|
(B9) |
|
(B10) |
把(B1) 式代入到(B6) 式至(B10) 式中，再根据(B5) 式就可以得到每一次的参考速度更正值.
当参考速度的更正值小于给定的值(本文选取10-6) 时，则终止迭代，并输出此时的参考速度.输出的参考速度可以使AWWE的精度在所选择的角度范围内达到全局最优.也就是我们要求的优化参考速度.
致谢感谢中海油研究总院地球物理总师李绪宣、开发研究院院长胡光义的指导；感谢中海油天津分公司、中海油服提供的VSP资料；感谢重大专项项目组及联合单位中国石油大学(北京) 提供的帮助.
| CaiX H, Liu Y, Wang J M, et al. 2015. Full-wavefield VSP reverse-time migration based on the adaptive optimal finite-difference scheme. Chinese J. Geophys.(in Chinese) , 58 (9) : 3317-3334. DOI:10.6038/cjg20150925 | |
| Chen H M, Zhou H, Lin H, et al. 2013. Application of perfectly matched layer for scalar arbitrarily wide-angle wave equations. Geophysics , 78 (1) : T29-T39. DOI:10.1190/geo2012-0062.1 | |
| Claerbout J F. 1976. Fundamentals of Geophysical Data Processing. Oxford:Scientific Publications. | |
| Guddati M N, Heidari A H. 2005. Migration with arbitrarily wide-angle wave equations. Geophysics , 70 (3) : S61-S70. DOI:10.1190/1.1925747 | |
| Guddati M N. 2006. Arbitrarily wide-angle wave equations for complex media. Computer Methods in Applied Mechanics and Engineering , 195 (1-3) : 65-93. DOI:10.1016/j.cma.2005.01.006 | |
| He B S, Zhang H X, Zhang J. 2008a. Prestack reverse-time depth migration of arbitrarily wide-angle wave equations. Acta Seismologica Sinica(in Chinese) , 30 (5) : 491-499. | |
| He B S, Zhang H X, Han L H, et al. 2008b. Prestack reverse-time depth migration of two acoustic wave equation and comparison between both. Oil Geophysical Prospecting (in Chinese) , 43 (6) : 661-684. | |
| He B S, Zhang H X, Zhang J J. 2010. Absorbing boundary conditions of arbitrarily wide-angle acoustic equations reverse-time migration. Journal of China Coal Society (in Chinese) , 35 (1) : 106-109. | |
| Heidari A H, Guddati M N. 2006. Highly accurate absorbing boundary conditions for wide-angle wave equations. Geophysics , 71 (3) : S85-S97. DOI:10.1190/1.2192914 | |
| Li W J, Qu S L, Wei X C, et al. 2012. Discussion about reverse-time depth migration technology of prestack elastic wavefield for offset VSP data. Chinese J. Geophys.(in Chinese) , 55 (1) : 238-251. DOI:10.6038/j.issn.0001-5733.2012.01.023 | |
| Ma Z T. 1983. A splitting-up method for solution of higher-order migration equation by finite-difference scheme. Acta Geophysica Sinica (in Chinese) , 26 (4) : 377-388. | |
| Sun Q F, Du Q Z. 2011. The frequency-space domain prestack depth migration with arbitrarily wide-angle wave equation. Oil Geophysical Prospecting (in Chinese) , 46 (6) : 890-896. | |
| Sun W B, Sun Z D. 2010. VSP reverse time migration based on the pseudo-spectral method and its applications. Chinese J. Geophys.(in Chinese) , 53 (9) : 2196-2203. DOI:10.3969/j.issn.0001-5733.2010.09.020 | |
| Zhou H, Chen H M, Li X X, et al. 2014. Methods for improving computational efficiency of arbitrarily wide-angle wave equation and suppressing evanescent wave noises. Chinese J. Geophys. (in Chinese) , 57 (4) : 1188-1198. DOI:10.6038/cjg20140416 | |
| 蔡晓慧, 刘洋, 王建民, 等. 2015. 基于自适应优化有限差分方法的全波VSP逆时偏移. 地球物理学报 , 58 (9) : 3317–3334. DOI:10.6038/cjg20150925 | |
| 何兵寿, 张会星, 张晶. 2008a. 任意广角波动方程叠前逆时深度偏移. 地震学报 , 30 (5) : 491–499. | |
| 何兵寿, 张会星, 韩令贺, 等. 2008b. 两种声波方程的叠前逆时深度偏移及比较. 石油地球物理勘探 , 43 (6) : 661–684. | |
| 何兵寿, 张会星, 张建军. 2010. 任意广角声波方程逆时偏移的吸收边界条件. 煤炭学报 , 35 (1) : 106–109. | |
| 李文杰, 曲寿利, 魏修成, 等. 2012. 非零偏VSP弹性波叠前逆时深度偏移技术探讨. 地球物理学报 , 55 (1) : 238–251. DOI:10.6038/j.issn.0001-5733.2012.01.023 | |
| 马在田. 1983. 高阶方程偏移的分裂算法. 地球物理学报 , 26 (4) : 377–388. | |
| 孙歧峰, 杜启振. 2011. 任意广角波动方程频率-空间域叠前深度偏移成像. 石油地球物理勘探 , 46 (6) : 890–896. | |
| 孙文博, 孙赞东. 2010. 基于伪谱法的VSP逆时偏移及其应用研究. 地球物理学报 , 53 (9) : 2196–2203. DOI:10.3969/j.issn.0001-5733.2010.09.020 | |
| 周辉, 陈汉明, 李绪宣, 等. 2014. 提高任意广角波动方程计算效率及压制倏逝波干扰方法. 地球物理学报 , 57 (4) : 1188–1198. DOI:10.6038/cjg20140416 | |