刘炜，王彦春，毕臣臣，徐仲博

(1. 成都理工大学 地球物理学博士后科研流动站，四川 成都 610059；2.中国地质大学(北京) 地球物理与信息技术学院，北京 100083)

0 引言

近年来，随着油气资源勘探开发程度的不断加深，勘探对象也日益复杂化，因此薄互层、高陡构造、微型构造、隐蔽性油气藏和岩性油气藏等复杂地质目标的精确识别也愈发重要。相比于常规地面地震观测方式，VSP观测方式一般采用地表激发、井中接收的地震数据采集方式，因此VSP资料具有波场信息丰富、分辨率和信噪比高以及环境噪声小等特点，有利于探明地下复杂地质构造以及提取储层弹性参数等信息[1]。逆时偏移方法是一种基于双程波波动方程的偏移成像方法，其不受地层倾角限制同时能够适应任意复杂速度区域，因此被认为是目前成像精度最高的地震资料偏移成像方法[2]。VSP地震数据和逆时偏移方法的有机结合，有利于精确识别地下复杂构造，因此研究VSP地震数据的逆时偏移方法具有一定的现实意义。

Whitemore于1983年首次提出了逆时偏移方法[3]，同年其他一些学者也对该方法进行了详细的探讨[4-5]。经过三十多年的发展，逆时偏移方法的理论体系逐渐完善，应用对象和范围也逐渐广泛，出现了叠后逆时偏移和叠前逆时偏移、声波逆时偏移[6]和弹性波逆时偏移[7]、各向同性介质逆时偏移[8]、各向异性介质逆时偏移[9]以及地面地震数据逆时偏移和VSP地震数据逆时偏移[10-12]等。逆时偏移的实现过程主要涉及到以下几个关键方面：波场延拓、边界条件、存储策略、成像条件以及低频噪声压制方法等。波场延拓是逆时偏移方法的核心，其实质是离散求解波动方程，目前有限差分法最为常用[13]。在利用有限差分法求解波动方程时，提高其数值模拟精度是需要重点考虑的内容，主要措施有增大差分阶数和寻求优化差分系数。根据频散关系求解有限差分法的差分系数时，通常有泰勒级数展开法和最优化方法两种方法，一般基于后者所得到的差分系数的数值模拟精度更高。Dablain基于泰勒级数展开法求取了高阶差分系数，有效地提高了有限差分法的数值模拟精度[14]。Liu和Sen利用泰勒级数法展开基于声波方程的频散关系分别获得了传统规则网格和交错网格有限差分法的差分系数[15-16]。Liu采用最小二乘法优化根据频散关系建立的目标函数分别求取了传统规则网格和交错网格有限差分法的优化差分系数，进一步提高了有限差分法的数值模拟精度[17-18]。在边界条件和存储策略方面，Clapp详细探讨了逆时偏移方法的存储需求，提出了边界存储策略，有效地减少了震源波场存储量[19]。在此基础上，王保利等提出了有效边界存储策略，进一步有效地减少了逆时偏移方法的存储需求[20]。段沛然等提出了基于优化算子边界存储策略的高效逆时偏移方法，其能够兼顾逆时偏移方法的存储量和计算时间[21]。在成像条件以及噪声压制方面，不少学者也进行了大量的研究。常用的逆时偏移方法成像条件有：反褶积成像条件[22]、零延迟互相关成像条件[23]、震源归一化零延迟互相关成像条件[24]等；常用的低频噪声压制方法有：拉普拉斯滤波法[25]和高阶拉普拉斯滤波法[26]等。

本文在二维变密度声波波动方程的基础上，基于优化交错网格有限差分法研究了VSP地震数据的高精度逆时偏移方法。优化交错网格有限差分法采用最小二乘法来确定其差分系数，因此相比于传统的交错网格有限差分法，当差分算子长度相等时，它能够更有效地压制数值频散，从而获得更高精度的地震波场数值模拟结果。为了进一步提高VSP数据逆时偏移方法的成像精度和实用性，采用PML吸收边界条件和有效边界存储策略在有效地压制地震波场边界反射的同时极大地减少震源波场的存储需求，采用震源归一化零延迟互相关成像条件进行成像计算，同时采用高阶拉普拉斯滤波法压制成像结果的低频成像噪声，进而获得高精度的成像结果。不同的模型测试结果表明：VSP地震数据逆时偏移方法可以实现对复杂地质构造的精确成像；由于VSP地震资料具有波场信息丰富、分辨率和信噪比高等特点，因此相比于常规地面地震数据的逆时偏移方法，VSP数据逆时偏移对井旁构造、高陡构造、微型构造等复杂地质构造进行成像更具有优势。

1 方法原理

逆时偏移方法的实现过程大体可以分为两步：首先进行震源波场的正向延拓和检波点波场的逆向延拓，然后采用成像条件对震源波场和检波点波场进行成像计算，进而得到最终的逆时偏移成像结果[27]。逆时偏移实现过程一般会涉及到以下几个关键性问题：波场延拓、边界条件、存储策略、成像条件[28]以及低频噪声压制方法等。

1.1 优化交错网格有限差分法

本文采用优化交错网格有限差分法离散求解变密度声波波动方程实现逆时偏移方法的波场延拓过程。二维变密度声波波动方程可以表示为[29]：

(1)

式中:ρ为介质密度；p为波场值；K=ρv2为体积模量，v为介质速度。

采用交错网格有限差分法对式(1)左右两边进行离散求解可得[30]：

(2)

式中：Δt为时间步长；h为空间网格步长，本文设定x和z方向的空间网格步长相等；M为差分算子长度，其值等于差分阶数的一半；am为差分系数。

将式(2)代入式(1)，根据链式法则并进行简化可得到二维变密度声波波动方程时间2阶、空间2M阶的交错网格差分格式[31]：

(3)

式中，r=vΔt/h；al和am为差分系数。

为了进一步提高波场延拓的模拟精度，采用基于最小二乘法的优化交错网格有限差分法来计算每个时刻的波场值。该方法主要是采用最小二乘法优化根据空间域频散关系构建得到的目标函数从而得到优化差分系数，其具体计算公式为[18]：

(n=2,3,…,M)

(5)

其中:

ψm(β)=2{sin[(m-0.5)β]-

2(m-0.5)sin(0.5β)}g(β)=β-2sin(0.5β)。(6)

式中：β=kh，k为波数；b为介于(0,π)之间的常数。

与传统基于泰勒级数展开法的交错网格有限差分法相比，基于最小二乘法的优化交错网格有限差分法可以进一步压制数值频散，提高地震波场模拟精度，尤其是在大波数范围内。传统交错网格有限差分法的差分系数计算公式为[16,31]：

(7)

式中，m=1,2,…,M。为了比较分析本文优化交错网格有限差分法的优越性，定义下列两个参数[31]：

(8)

其中:

(9)

式中：δ(β,θ)为二维声波波动方程数值模拟的相速度频散；vFD为有限差分数值求解的相速度；θ为平面波的传播方向；ε(β,θ)为二维正演模拟时一个网格间距内地震波传播时间的相对误差。由式(8)可知，当δ越接近于1或ε越接近于0时，二维地震波场数值模拟的数值频散就越小，即模拟精度越高。

为了比较分析传统和优化交错网格有限差分法的数值频散特征，设计一个均匀模型：介质速度为 1 500 m/s，密度为2.0 g/cm3，模型大小为4 000 m×4 000 m，空间网格步长为20 m，时间步长为1 ms。图1分别展示了不同差分算子长度时传统和优化交错网格有限差分法的数值频散误差随波数的变化规律。从图中可以看出，随着差分算子长度增加，基于优化和传统交错网格有限差分法的地震波场数值模拟的频散误差均会减小，即模拟精度均会提高；在差分算子长度相等时，优化交错网格差分法的模拟精度优于传统交错网格差分法。图2分别展示了相同差分算子长度时传统和优化交错网格有限差分法的频散误差随地震波传播方向的变化规律。从图中可以看出，对于两种交错网格差分法，地震波场的模拟精度均与地震波的传播方向有关；在地震波传播方向相同时，优化交错网格差分法的模拟精度优于传统交错网格差分法。设定震源子波是主频为15 Hz的Ricker子波且位于模型正中心，地震记录总时长为3 s，对该模型进行正演模拟进而对比分析两种交错网格差分法的模拟效果。图3为不同差分算子长度时传统和优化交错网格差分法在不同时刻的波场快照。从图中可以看出，随着差分算子长度增大，交错网格差分法的数值频散降低，模拟精度提高；在相同的差分算子长度条件下，优化交错网格差分法具有更高的数值模拟精度。

a—传统交错网格有限差分法;b—优化交错网格有限差分法a—conventional staggered-grid finite-difference;b—optimal staggered-grid finite-difference图1 不同差分算子长度时交错网格有限差分法的频散误差随波数的变化规律Fig.1 Variations of dispersion error with wavenumber and different operator lengths by using different staggered-grid finite-difference methods

a—传统交错网格有限差分法;b—优化交错网格有限差分法a—conventional staggered-grid finite-difference;b—optimal staggered-grid finite-difference图2 相同差分算子长度时(M=5)交错网格有限差分法的频散误差随传播方向的变化规律Fig.2 Variations of dispersion error with different propagation angles by using different staggered-grid finite-difference methods (M=5)

a—传统交错网格有限差分法，M=5;b—优化交错网格有限差分法，M=5;c—传统交错网格有限差分法，M=10;d—优化交错网格有限差分法，M=10；从左至右依次为1 s和2.5 sa—conventional staggered-grid finite-difference for M=5;b—optimal staggered-grid finite-difference for M=5;c—conventional staggered-grid finite-difference for M=10;d—optimal staggered-grid finite-difference for M=10;from left to right time at 1s and 2.5s图3 不同差分算子长度时交错网格有限差分法在不同时刻的波场快照Fig.3 Snapshots by using different staggered-grid finite-difference methods and different operator lengths

算法稳定性也是衡量有限差分数值模拟方法优劣的重要因素之一。交错网格有限差分法的稳定性条件可以表示为[16,31]：

(10)

式中s为稳定性因子，其值越接近于1表示有限差分算法稳定性越好。

利用上述均匀模型对本文交错网格有限差分法进行稳定性分析。图4给出了基于泰勒级数展开法的传统交错网格有限差分法以及基于最小二乘法的优化交错网格有限差分法的稳定性因子曲线。由图可知，优化交错网格有限差分法的稳定性略差于传统交错网格有限差分法。通过以上分析可知，优化交错网格有限差分法在满足稳定性的条件下更加有利于压制数值频散，因此本文在后续的逆时偏移方法研究过程中主要采用时间2阶、空间16阶的优化交错网格有限差分法进行波场延拓。

图4 不同交错网格有限差分法的稳定性因子曲线Fig.4 Curves of stability factor of different staggered-grid finite-difference methods

1.2 PML边界条件及有效边界存储策略

在地震波场的数值模拟过程中，由于计算机的计算区域有限，在模型四周存在的人工截断边界会造成强烈的边界反射从而影响数值模拟结果，因此在地震波场的数值模拟过程中需采用有效的边界条件来压制边界反射，从而提高地震波场数值模拟结果的准确性。PML吸收边界条件可以吸收任意方向、任意频率的波，因此目前被广泛应用于地震波场的数值模拟过程[32-33]。

二维变密度声波波动方程的PML吸收边界条件的控制方程可以表示为：

(11)

式中，A为吸收衰减系数矩阵，其值从PML吸收边界的内边界开始由内到外依次增大，具体取值由吸收衰减因子决定。式(11)的交错网格离散格式为：

(12)

PML吸收边界条件的衰减因子分布情况如图5所示。图中区域E为实际模型区域，吸收衰减因子为零；区域A、B、C、D、F、G、H、I为吸收衰减区域，吸收衰减因子不为零。吸收衰减因子的种类繁多，本文采用余弦型衰减因子[34]，其具体表达式为：

(13)

式中：βx和βz分别表示沿x和z方向的余弦型衰减因子；B为衰减幅度因子；Lx和Lz分别表示x和z方

图5 PML吸收边界条件简易示意Fig.5 Simple sketch of PML absorbing boundary condition

向完全匹配层的总网格数；lx和lz分别表示x和z方向距离完全匹配层外边界的网格数。

PML吸收边界条件能够有效地压制由人工截断边界造成的边界反射，从而极大地削弱边界反射对地震波场数值模拟结果中有效信息的影响，因此能够有效地满足逆时偏移成像的精度要求。然而，逆时偏移方法要求事先已知每个时刻的震源波场，存储所有时刻的震源波场对计算机的存储要求很高，一般计算机无法满足。因此本文在PML吸收边界条件的基础上引进有效边界存储策略，以求在保证波场延拓精度的条件下极大地减少逆时偏移震源波场的存储需求。由图5可知，在波场延拓的过程中，实际的正演模型包含内部有效模型区域和外部吸收衰减边界区域，即如图6所示，区域A(深灰色)为内部有效模型区域，区域B(白色)和区域C(浅灰色)为PML吸收边界区域。有效边界存储策略只需存储每个时刻区域C的震源波场值以及最后两个时刻所有区域的震源波场值就可以精确地逆时重构出所有时刻的震源波场，从而有效地降低逆时偏移的震源波场存储需求。区域C的厚度与差分算子长度有关，在本文中波场延拓采用变密度声波波动方程，因此区域C的网格点数等于差分算子长度的两倍，即等于差分阶数。

图6 有效边界存储策略简易示意Fig.6 Simple sketch of effective boundary storage strategy

为了验证上述PML吸收边界条件以及基于有效边界存储策略重构震源波场的正确性和有效性，设计如图7所示的多层层状模型。该模型大小为4 000 m×4 000 m，空间网格步长为10 m，时间步长为1 ms，地震记录总时长为3 s，PML吸收边界的网格数为50，震源子波是主频为15 Hz的Ricker子波且位于(2 000 m，1 800 m)处。图8分别展示了该模型在不同时刻的正传波场快照、应用有效边界存储策略时所重构的相对应时刻的震源重构波场快照以及它们二者之间的差异剖面。从图中可以看出，在地震波场正向延拓过程中，当地震波传播至模型边界时，地震波会继续传播到PML吸收边界内，然后PML吸收边界对这些地震波进行吸收衰减，因此由人工截断边界造成的边界反射得到了极大地削弱，说明PML吸收边界条件可以有效地压制边界反射，降低其对有效信号的影响。对比分析图8a、图8b和图8c可以发现，在同一时刻震源的正传波场和重构波场的波形基本一致，二者之间的数值差异也很小，证明了采用有效边界存储策略在减小震源波场存储需求的同时能够根据所存储的部分震源波场信息精确地重构出震源波场。综上可知，在逆时偏移方法中，综合应用PML吸收边界条件以及有效边界存储策略，可以在有效压制边界反射的同时极大地降低震源波场的存储需求。

图7 多层层状模型Fig.7 A multilayer model

a—正传波场;b—重构波场;c—正传波场和重构波场之间的差异;从左至右依次为0.4、0.8、1.2 s时刻a—forward wavefields;b—reconstructed wavefields;c—difference between the forward and reconstructed wavefields;from left to right time at 0.4,0.8,1.2 s图8 多层层状模型在不同时刻的波场快照Fig.8 Snapshots at different time for the multilayer model

1.3 成像条件

成像条件是影响逆时偏移成像精度的关键因素之一。目前，逆时偏移方法的成像条件种类繁多，如激发时间成像条件、最大振幅到达时成像条件、反褶积成像条件和零延迟互相关成像条件等。其中，零延迟互相关成像条件由于既能充分利用全部地震波场信息又能在一定程度上压制成像噪声，因此被广泛应用于逆时偏移方法中。其数学表达式为[23]：

(14)

式中：I(x,z)为成像值；S(x,z,t)为震源波场；R(x,z,t)为检波点波场；Tmax为地震记录总长度。

由式(14)可知，零延迟互相关成像条件通过对震源波场和检波点波场进行互相关计算从而获得最终的逆时偏移成像结果，但是无法削弱震源对成像结果的影响，同时最终成像结果也无法直接反映地下地层的反射系数。在此基础上，有学者提出了震源归一化零延迟互相关成像条件，能够有效地弥补零延迟互相关成像条件的缺陷，进一步提高逆时偏移的成像精度。震源归一化零延迟互相关成像条件的具体数学表达式为[24]：

(15)

对比式(14)和式(15)可知，震源归一化零延迟互相关成像条件相当于在零延迟互相关成像条件的基础上除以震源波场的互相关结果，因此能够有效地削弱震源对成像结果的影响，同时所得成像值可以直接对应于地下地层的反射系数。

1.4 噪声压制

由于逆时偏移方法本身的原因，在其成像结果中不可避免地会产生低频成像噪声，为了提高成像精度，需采用有效方法对其进行压制。拉普拉斯滤波法是一种常用的噪声压制方法，相当于将一个滤波器直接作用于逆时偏移成像数据，具有简单易实施、噪声压制效果较好等特点。但是经过一些学者的研究证明，常规拉普拉斯滤波法无法完全地压制逆时偏移的成像噪声，在去噪后的剖面上仍然残留部分低频噪声。在前人的研究基础上，本文采用高阶拉普拉斯滤波法来压制逆时偏移的低频成像噪声，其表达式为[26]：

式中：N=4,6,…,2n，即为大于2的偶数。该公式的简易离散形式可以表示为：

式中：f(x,y)表示二维数据在(x,y)的数值。式(17)的上角标并非表示指数，而是表示拉普拉斯算子的阶数。

利用傅里叶变换将式(16)转换到波数域可得：

(18)

式中：ω为角频率，θ为地震波入射角。从式(18)可以看出，拉普拉斯类滤波法相当于对逆时偏移成像数据进行角度域衰减，其衰减程度与拉普拉斯算子的阶数有关，不同阶数的拉普拉斯滤波法仅衰减系数存在差异，分析可知高阶拉普拉斯滤波法对大角度的逆时偏移成像噪声有更好的压制效果，但是阶数过大也会加大损害逆时偏移结果中的有效信号，因此本文令N值等于4。

2 模型测试

为了验证本文VSP地震数据逆时偏移方法的有效性，采用两个模型进行测试，分别为二维SEG/EAGE盐丘模型和Marmousi模型。

2.1 二维SEG/EAGE盐丘模型

二维SEG/EAGE盐丘模型的速度和密度如图9所示，其模型大小为3 380 m×2 100 m，空间网格步长为10 m，时间步长为1 ms，震源为15 Hz的Ricker子波，地震记录总长度为4 s。对于该模型，采用常规地面地震观测方式和VSP观测方式进行逆时偏移方法测试，分别对比分析了它们的逆时偏移成像结果。常规地面地震观测系统的参数为：震源均匀地分布在地表，炮间距为30 m，共113炮；检波器同样均匀地分布在地表，道间距为10 m，共338道。VSP观测系统参数为：震源分布方式与常规地面地震观测系统一致，但是其检波器则均匀地分布在井口坐标为(0 m，0 m)和(3 380 m，0 m)的两口井中，道间距为10 m，两口井共420道。针对该SEG/EAGE盐丘模型，主要分析讨论了零延迟互相关成像条件与震源归一化零延迟互相关成像条件的成像效果、高阶拉普拉斯滤波法的低频成像噪声压制效果(包括其对成像结果振幅和相位的影响)以及常规地面地震数据逆时偏移和VSP数据逆时偏移的差异。

图9 二维SEG/EAGE盐丘模型Fig.9 2D SEG/EAGE salt model

图10显示了分别应用零延迟互相关成像条件和震源归一化零延迟互相关成像条件时该二维SEG/EAGE盐丘模型第57炮的VSP数据逆时偏移成像结果。从图中可以看出，在基于零延迟互相关成像条件的逆时偏移结果中,震源对其成像结果的影响较大，而在基于震源归一化零延迟互相关成像条件的逆时偏移结果中,震源对其成像结果的影响得到了极大的削弱，说明震源归一化零延迟互相关成像条件能够有效地削弱震源对逆时偏移成像结果的影响，从而提高其成像结果的精度。图11给出了应用震源归一化零延迟互相关成像条件时采用高阶拉普拉斯滤波法进行噪声压制前后该二维SEG/EAGE盐丘模型的地面地震数据和VSP数据逆时偏移结果。从图中可以看出，在噪声压制前，常规地面地震数据和VSP地震数据的逆时偏移结果中均存在明显的低频成像噪声，采用高阶拉普拉滤波法进行噪声压制后，二者的成像剖面均变得更为清晰，即成像精度提高，说明高阶拉普拉滤波法能够有效地压制逆时偏移结果的低频成像噪声。同时，相比于常规的地面地震数据逆时偏移结果，VSP数据的逆时偏移结果能够更加精确地反映地下介质的构造形态和构造特征，尤其是高陡构造部位以及井旁部位，说明VSP数据的逆时偏移方法更有利于精确地识别地下复杂地质构造。另外，由拉普拉斯算子的原理可知，拉普拉斯类滤波法实质是对成像结果的数据进行角度域滤波处理，其会对成像结果的振幅和相位产生一定的破坏。图12给出了应用高阶拉普拉斯滤波法前后该二维SEG/EAGE盐丘模型VSP逆时偏移结果的幅值谱和相位谱，其横纵坐标均表示空间频率。从图中可以看出，高阶拉普拉斯滤波法可以压制逆时偏移结果的成像噪声，噪声的振幅得到了有效的削弱，但是同时也可以发现滤波前后逆时偏移成像结果中有效信号的幅值谱和相位谱发生了改变，说明高阶拉普拉滤波法会在一定程度上损害逆时偏移结果的振幅信息和相位信息。

a—零延迟互相关成像条件;b—震源归一化零延迟互相关成像条件a—cross correlation imaging condition;b—normalized cross correlation imaging condition of sources图10 二维SEG/EAGE盐丘模型的第57炮VSP数据逆时偏移结果Fig.10 RTM results of the 57th VSP gather for the 2D SEG/EAGE salt model

a—噪声压制前幅度谱;b—噪声压制后幅度谱;c—噪声压制前相位谱;d—噪声压制后相位谱a—amplitude spectrum before noise suppression;b—amplitude spectrum after noise suppression;c—phase spectrum before noise suppression;d—phase spectrum after noise suppression图12 二维SEG/EAGE盐丘模型VSP逆时偏移结果的幅值谱和相位谱Fig.12 Amplitude spectrum and phase spectrum of VSP RTM results for the 2D SEG/EAGE salt model

2.2 Marmousi模型

为了进一步验证本文VSP逆时偏移方法对于复杂速度模型的有效性，采用Marmousi模型对其进行测试。图13给出了该Marmousi模型的速度和密度参数，其模型大小为7 000 m×7 000 m，空间网格步长为20 m，时间步长为1 ms，震源子波为15 Hz的Ricker子波，地震记录总长度为6 s。对于该Marmousi模型，本文同样对比分析了其常规地面地震数据和VSP地震数据的逆时偏移结果。两种观测系统的震源均以60 m为间距均匀分布在地表，共117炮，常规地面地震观测系统的检波器以20 m为间距均匀分布在地表，共350道，而VSP观测系统的检波器则以20 m为间距均匀地分布在井口坐标为(0 m，0 m)和(7 000 m，0 m)的两口井中，两口井共700道。针对该Marmousi模型，首先分析讨论了两种不同的互相关成像条件的成像效果、高阶拉普拉斯滤波法的低频成像噪声压制效果、常规地面地震数据逆时偏移和VSP数据逆时偏移的差异，另外还进一步讨论了逆时偏移方法对偏移速度的敏感性。

图14展示了分别应用零延迟互相关成像条件和震源归一化零延迟互相关成像条件时该Marmousi模型第59炮的VSP数据逆时偏移成像结果。从图中可以看出，与二维SEG/EAGE盐丘模型的单炮逆时偏移结果类似，震源归一化零延迟互相关成像条件极大地削弱了震源对VSP数据逆时偏移成像结果的影响，从而提高了VSP数据单炮逆时偏移成像结果的精度。

图13 Marmousi模型Fig.13 Marmousi model

图15显示了应用震源归一化零延迟互相关成像条件时采用高阶拉普拉斯滤波法进行噪声压制前后该Marmousi模型的地面地震数据和VSP数据逆时偏移结果。从图中可以看出，与二维SEG/EAGE盐丘模型的最终逆时偏移结果类似，无论是常规地面地震数据的逆时偏移结果还是VSP数据的逆时偏移结果，高阶拉普拉滤波法均有效地压制了逆时偏移结果中的低频成像噪声，提高了逆时偏移成像结果的精度。与此同时，VSP数据的逆时偏移结果比常规地面地震数据的逆时偏移结果更加精确，能够更加有效地预测地下复杂地质构造。综上可知，VSP数据的逆时偏移方法更加有利于识别和预测地下复杂地质目标。

a—零延迟互相关成像条件;b—震源归一化零延迟互相关成像条件a—cross correlation imaging condition;b—normalized cross correlation imaging condition of sources图14 Marmousi模型的第59炮VSP数据逆时偏移结果Fig.14 RTM results of the 59th VSP gather for the Marmousi model

a—地面数据低频噪声压制前;b—VSP数据低频噪声压制前;c—地面数据低频噪声压制后;d—VSP数据低频噪声压制后a—surface data before low frequency noise suppression;b—VSP data before low frequency noise suppression;c—surface data after low frequency noise suppression;d—VSP data after low frequency noise suppression图15 Marmousi模型的常规地面地震数据和VSP数据的逆时偏移结果Fig.15 RTM results of conventional surface seismic data and VSP data for the Marmousi model

为了分析本文逆时偏移方法对速度的敏感性，将该Marmousi模型的速度模型和密度模型进行平滑处理，其结果如图16a所示，然后利用该平滑后的模型进行逆时偏移，其结果如图16b所示。从图中可以看出，当偏移速度模型不准确时，常规地面地震数据逆时偏移结果以及VSP数据逆时偏移结果的精确度均会降低，尤其是在速度变化剧烈的部位，说明偏移速度是影响逆时偏移结果精度的重要因素之一。同时可以发现，此时VSP数据的逆时偏移结果仍然比地面地震数据的逆时偏移结果更为精确。

综上说明，即使偏移速度出现细微误差，在相同的条件下，VSP数据逆时偏移也能更为精确地识别地下地质构造。

a—速度;b—密度;c—地面地震数据逆时偏移结果;d—VSP数据逆时偏移结果a—velocity;b—density;c—RTM results of conventional surface seismic data;d—RTM results of VSP data图16 平滑后的Marmousi模型及其逆时偏移结果Fig.16 The smoothed Marmousi model and its corresponding RTM results

3 结论与认识

在前人研究的基础上，研究了基于优化交错网格有限差分法的VSP数据逆时偏移方法。针对逆时偏移的不同方面，采用了不同的策略和措施。对于波场延拓，采用基于最小二乘法的优化交错网格有限差分法进行波场的数值模拟，有效提高了数值模拟精度；采用PML吸收边界条件和有效边界存储策略相结合的方式，在极大压制边界反射的同时有效地降低逆时偏移震源波场的存储需求，转而采用波场重构的方式精确地重构震源波场。在成像的过程中，应用震源归一化零延迟互相关成像条件进行VSP数据的成像过程，有效地削弱震源对逆时偏移成像结果的影响。最后，采用高阶拉普拉斯滤波法压制逆时偏移成像结果的低频噪声，提高成像结果的精度。通过不同的模型测试，得到了以下几点结论和认识：

1)相比于传统的基于泰勒级数展开法的交错网格有限差分方法，基于最小二乘法的优化交错网格有限差分法能够有效地提高地震波场的数值模拟精度，尤其是在大波数范围内；交错网格有限差分法的数值模拟精度随着差分算子长度的增大而提高，随着地震波传播方向的变化而发生改变。在相同的条件下，优化交错网格有限差分法的模拟效果优于传统交错网格有限差分法。

2)PML吸收边界条件能够有效压制由于计算空间有限造成的边界反射；在此基础上，应用有效边界存储策略可以有效降低逆时偏移方法震源波场的存储需求，其重构震源波场与正传震源波场基本一致，满足高精度逆时偏移的精度要求。

3)震源归一化零延迟互相关成像条件能够有效地削弱震源对逆时偏移结果的影响，同时高阶拉普拉斯滤波法能够有效压制低频成像噪声，从而提高逆时偏移成像结果的精度，但是高阶拉普拉斯滤波法会在一定程度上破坏逆时偏移结果的振幅信息和相位信息。

4)相比于传统的地面地震数据逆时偏移，VSP数据的逆时偏移能够更加精确地识别地下复杂地质构造，尤其是井旁构造、高陡构造、微型构造以及速度变化剧烈构造等。