基于正态分布权函数的VSP-CDP叠加成像

2020-03-02杨飞龙李辉峰

石油地球物理勘探 2020年1期

杨飞龙李辉峰孙辉张雪罗浩赵驰

(①西安石油大学地球科学与工程学院，陕西西安710065； ②陕西省油气成藏地质学重点实验室，陕西西安710065；③山东科技大学地球科学与工程学院山东省沉积成矿作用与沉积矿产重点实验室，山东青岛 266590；④西南交通大学地球科学与环境工程学院，四川成都 611756)

0 引言

随着油气勘探的深入，面对规模小、结构复杂、物性变化剧烈的地层圈闭和岩性圈闭，对地震勘探分辨率的要求越来越高。垂直地震剖面法勘探(VSP)是一种井中地震观测技术，主要研究地球介质内部波场特征、地震参数以及井周围地质构造。与地面地震相比，VSP资料的分辨率更高，可以得到更精确的时深转换结果及速度模型，并为零相位子波分析提供支持。

前人针对VSP反射波资料使用基于射线理论的CDP叠加方法进行成像。Wyatt[1]首次提出VSP-CDP的概念，假设地层为常速水平层，当面对倾斜地层及复杂构造时，成像效果较差。Robert等[2]提出共中心深度点(CMD)叠加成像方法，基于水平反射界面假设，将共炮点道集数据转换成CMD道集数据进行叠加成像，利用多道叠加提高了信噪比，但只适合反射界面倾角较小的情况。Smalley等[3]提出共侧向点(CLP)叠加成像方法，假设在各向同性介质中共炮点道集满足双曲线时距方程，不适合横向速度变化剧烈的地层。蒋振武等[4]研究了针对斜层的井间地震共深度点(DLCDP)成像算法。严又生等[5]阐述了非均匀介质井间地震VSP-CDP反射波叠加成像算法。石星等[6]应用共反射面元叠加成像方法对复杂构造成像，信噪比较常规叠加剖面更高。邓金华等[7]提出了一种改进的共反射面元叠加方法，有效地避免了零井源距剖面中的同相轴相交现象。杨飞龙等[8]利用高斯束方法对VSP数据成像，提出了VSP高斯射线束法叠加成像算法。基于常规射线追踪理论的叠加成像方法主要解决VSP简单波场成像问题，面对复杂构造时成像精度较低。

人们将基于偏移理论的成像方法引入VSP波场成像。Authur等[9]利用层析与相位屏法波动方程偏移联合算法考察井中地震反射波成像。王华忠等[10]利用三维声波方程实现了VSP 数据叠前深度偏移，并给出对应方程及其差分格式。宋炜等[11]使用Kirchhoff积分法进行井间地震反射波偏移成像。Liu等[12]以Kirchhoff偏移为基础，使用波动路径叠前偏移成像方法对井中地震数据成像。Fei 等[13]利用相移加插值(Phase-Shift Plus Interpolation，PSPI)法对二维、三维 VSP 数据偏移成像。刘诗竹[14]详细研究了VSP数据的数值模拟和偏移成像方法。王维红等[15]对Walkaway VSP逆时偏移进行了界面集成，推动了VSP 资料处理的实用化进程。蔡晓慧等[16]提出了基于自适应优化有限差分法的VSP 逆时偏移算法，在保证算法精度的前提下，提升了计算效率。杨继东等[17]在高斯束偏移的基础上，通过修改射线束传播算子，提出了菲涅耳束偏移方法。杨飞龙[18]将高斯束偏移方法扩展到井中地震勘探，实现了井中地震高斯束叠前深度偏移成像方法。王冲等[19]利用Walkaway-VSP技术处理吐哈油田三塘湖盆地VSP资料，详细刻画了火山岩的形态和展布范围。俞岱等[20]使用高斯束偏移方法试算井中地震模型及实际资料，通过对比成像结果与地面地震剖面验证了方法的稳定性。基于偏移理论的成像技术受观测系统限制，导致成像剖面边缘的“画弧”现象严重。

在射线类偏移成像中，权函数对成像精度有一定影响。Liu等[21]研究了束偏移中的权函数，提出了适合Kirchhoff束偏移的叠加权函数。Sun等[22-23]在文献[21]的基础上改进权函数计算方法，提出了基于压缩感知的Kirchhoff束偏移方法。若使用上述权函数对VSP数据进行VSP-CDP转换，会出现波场能量失真，不能实现保幅叠加成像。

为此，本文使用动态射线追踪方法模拟VSP复杂构造地震波场，通过交互对比正演波场与实际波场特征获取准确的构造模型及速度场，最后使用正态分布权函数对VSP波场进行VSP-CDP转换。前人在VSP-CDP转换时仅根据地震波的传播特征将共炮点道集数据中的每道数据归位到相应的反射点处，本文使用正态分布权函数归位VSP波场的同时，将深度—时间域的每一个采样点转换成反射点井源距—深度域的多个样点，使反射点的覆盖次数均匀，提高了VSP成像精度。

1 理论方法

1.1 动态射线追踪

在直角坐标系下，波动方程可以表示为

(1)

式中：u为波场；v为传播速度；t为传播时间；x和z为笛卡尔坐标。

射线追踪将波场分解为射线场，在已知射线参数的条件下，可选任一射线并建立相应的射线中心坐标系(s,n)(图1)。

图1 中心坐标系下射线传播示意图

对于任意一条射线，(s,n)是中心射线Ω上某一点的坐标，s代表地震波传播的距离，n代表相邻射线与Ω的垂直距离。坐标系的基矢量分别为与Ω相切的单位切向量t和与Ω垂直并指向Ω内侧的单位法向量n

把式(1)变换到射线中心坐标系(s,n)，有

(2)

由式(2)可得到集中于射线中心邻近的解。抛物型波动方程具有下列形式的时间调和解

(3)

式中ω为角频率。式(3)表示沿中心射线进行积分。

对式(3)进行高频近似并忽略ω的高阶项可得u(s,n,t)的近似表达式

(4)

其中

(5)

Γ,s+vΓ2+v-2v,nn=0

(6)

(7)

式中：Γ=Γ(s)为位置的复值函数；A为相邻射线与中心射线的权函数关系；Γ,s、A,s分别为Γ、A对s的一阶偏导数；v,nn为v对n的二阶偏导数。

式(6)为动态射线追踪方程，式(7)为传输方程。由式(6)和式(7)解出Γ(s)和A(s)后可求出式(5)的W(s,v)，从而求得抛物型波动方程的解。

引入一新的复值函数q(s)，设

(8)

式中q,s为q(s)对s的一阶偏导数。将式(8)代入式(6)，得到关于q的二阶线性微分方程

vq,ss-v,sq,s+v,nnq=0

(9)

式中q,ss为q(s)对s的二阶偏导数。再令q,s=vp，则式(9)可写为

v(v,sq+vp,s)-v,s(vp)+v,nnq=0

(10)

式中v,s、p,s分别为v、p对s的一阶偏导数。进一步化简，得到一阶微分方程组

(11)

将式(11)的第1式代入式(8)得

Γ=pq-1

(12)

将其变换到笛卡尔坐标系，得

(13)

(14)

式中：τ(s)为地震波旅行时；K(s)为地震波的波前曲率；Lh为地震波传播到检波点处的有效半宽度，与角频率及动态射线参数有关。则式(13)变为

(15)

任意位置(s，t)的地震波场值为

(16)

式中：φ为射线入射的角度；uφ(s,n)为在(s,n)处的地震波场值。

1.2 正态分布权函数

常规的VSP-CDP叠加成像中，仅将地震反射波场转换至相应地层的一个反射点处，难以满足炮数少、构造复杂情形下的成像精度。本文从动态射线追踪有效邻域波场近似理论出发，在合成地震记录时设检波点处的波场为Lh范围内所有射线的能量加权，因此进行VSP-CDP转换时考虑有效范围内所有射线能量的高斯加权，并根据权函数将检波点的波场分解至Lh范围内所有反射点处，大大提高了覆盖次数。Liu等[21]认为：射线束能量由中心射线往两边衰减，距中心射线越近，对应的振幅值越大；通过计算地震波走时发现，射线束内网格节点走时等相关信息是由泰勒近似求取的，距中心射线越近的点，相对误差越小，精度越高，对应的权重越大。因此，Liu等[21]提出使用与距离相关的余弦平方窗函数作为成像公式中权函数主体，即

(17)

式中ns为相邻射线与中心射线之间的距离。前人的研究结果以及动态射线有效邻域波场近似理论表明，距检波点越近的射线对能量的贡献越大(图2)。在VSP-CDP转换时为了定量描述每条射线对检波点能量的贡献，假设检波点处的波场值为1，则在Lh范围内相邻射线能量的加权等于1，并且距检波点越近的射线对能量贡献越大，符合正态概率密度函数分布。正态分布的概率密度函数与x轴围成的面积为1，即将Lh范围内的射线能量加权作为该检波点的能量。根据正态概率密度函数的定义

(18)

式中μ、σ(σ>0)为两个常数。σ越小，正态曲线越陡峭；σ越大，正态曲线越平坦。

参考正态曲线的形态及物理意义，可以将相邻射线与中心射线权函数关系A近似表示为

(19)

当相邻射线与中心射线重合时(ns=0)，由单一中心射线的波场能量表征该接收道的波场(常规射线追踪方法)；当nsLh时，相邻射线的波场能量对该接收道的波场无贡献。

图2 中心射线与相邻射线权函数关系示意图

利用正态分布权函数将一道地震记录进行VSP-CDP转换至Lh范围内有效反射点，那么有效反射点处能量的加权应为对应的检波点波场能量。设某一道地震波场值为1，且该地震道的Lh为100m，利用权函数对该道地震波场进行VSP-CDP转换，得到ns变化时反射点处波场特征示意图(图3)。可见，随着相邻射线与中心射线距离增大，地震波场值逐渐减小。若使用文献[21]的权函数进行VSP-CDP转换，则会出现有效射线的波场加权大于原始地震波场(图3a)；本文提出的正态分布权函数使有效射线的波场加权等于原始地震波场(图3b)，因此可保幅处理地震资料。

图3 ns变化时反射点处波场特征示意图

1.3 VSP-CDP转换

1.3.1 常规VSP-CDP转换原理

常速水平层状介质VSP反射波时距曲线可以表示为(图4)

(20)

式中：V为平均速度；x0为井源距；zR为检波点深度；t为反射波旅行时。

VSP-CDP转换的实质是把共炮点道集(或共检波点道集)数据的每一个采样点从深度—时间域变换到井源距—反射点时间(深度)域，形成共反射点道集数据。根据几何关系(图4)可知

(21)

式中(x1，z1)为反射点坐标。式(21)给出了常速水平地层情况下样点从(z，t)域到(x，z)域的转换公式。

图5为单道数据的VSP-CDP转换示意图。由图可见，根据上述VSP-CDP转换方法可以计算到达检波点R的两个反射界面的反射点坐标(x1,z1)和(x2,z2)，分别对应T1和T2时刻，VSP-CDP转换将共炮点道集记录中T1、T2时刻之间的地震记录分别归位到反射点(x1,z1)、(x2,z2)处(若要转换成井源距—反射点深度域，需要对T1时刻之前的地震记录进行时深转换，再进行反射点归位)，这样就形成了两个反射界面情况下VSP单道数据的VSP-CDP转换。同理，可完成单道数据所有反射界面及整个VSP多道数据的VSP-CDP转换。

图4 常速水平地层VSP-CDP转换几何示意图

图5 单道数据的VSP-CDP转换示意图

1.3.2 基于正态分布权函数的VSP-CDP转换原理

在进行VSP-CDP转换时，利用动态射线追踪正演模拟计算反射点位置及检波点有效邻域范围内射线对能量贡献的正态分布权函数，并采用该权函数对波场分离后的VSP野外数据进行VSP-CDP转换，将共炮点道集数据的一个样点分解至共反射点道集数据的多个样点，即将检波点R的波场根据相邻射线与中心射线正态分布权函数关系分解至Lh范围内的所有射线对应的反射点处(图6)。最后按照一定的间隔划分面元，并对同一面元内的样点进行叠加。对所有样点进行上述操作，就得到最终的叠加成像剖面。与常规的VSP-CDP叠加成像方法相比，基于正态分布权函数的VSP-CDP方法有效增加了相同叠加面元内的反射点个数，进行VSP-CDP转换后，有效扩大了反射点的照明范围，大大改善了VSP勘探的横向分辨率。

图6 基于正态分布权函数的VSP-CDP转换原理

2 模型试算

使用由表1参数表征的观测系统对地质模型(图7，模型参数如表2所示)进行VSP波场正演模拟，得到切除初至后的原始地震波场(图8)。根据地质及测井资料建立初始速度场，并在同一观测系统下进行VSP动态射线追踪正演数值模拟(图9)。调整速度及地质模型的过程是一个迭代过程，通过交互地质建模软件微调地质模型，同时根据地震反射同相轴斜率修改地层速度，直至正演波场与原始波场吻合为止，图10为第2炮数据原始波场与正演波场。使用基于正态分布权函数的VSP-CDP叠加成像方法进行反射P波和反射SV波成像，并选择10m×10m叠加面元得到叠加剖面(图11)。可见，基于正态分布权函数的VSP-CDP叠加成像方法对地下复杂构造及微小构造的成像精度较高。

表1 野外采集观测系统参数

图7 地质模型

表2 地质模型初始参数

图8 切除初至后的原始地震波场

图9 VSP动态射线追踪正演波场

图10 第2炮数据原始波场(黑色)与正演波场(红色)

图11 基于正态分布权函数的VSP-CDP叠加成像剖面

3 方法对比

为了验证基于正态分布权函数的VSP-CDP叠加成像方法的效果，对地质模型(图7)的局部区域使用常规射线追踪的VSP-CDP叠加、基于正态分布权函数的VSP-CDP叠加方法成像，选择10m×10m的横向叠加面元得到叠加剖面(图12)。可见，常规VSP-CDP叠加成像方法在靠近井旁构造区域及断层面以下区域出现空道(图12a)，基于正态分布权函数的VSP-CDP叠加成像方法没有出现空道(图12b)，表明前者的横向分辨率低于后者。

图12 不同成像方法效果对比

4 实际资料测试

为了验证基于正态分布权函数的VSP-CDP叠加成像方法的效果，对银额盆地哈日凹陷东部洼陷区Y5井非零井源距VSP资料进行成像。该VSP资料包括1炮零井源距数据、1炮非零井源距数据，图13为Y5井非零井源距VSP三分量数据，图14为Y5井非零井源距VSP三分量合成数据。根据研究区三维地面地震解释结果建立实际地质模型(图15)，通过零井源距VSP资料获取井旁地层速度，然后使用VSP动态射线追踪方法正演模拟，并将正演记录与VSP上行反射纵波波场交互对比、修改速度模型，直至两者吻合为止(图16)。分别使用常规VSP-CDP叠加方法和基于正态分布权函数的VSP-CDP叠加方法对Y5井非零井源距VSP资料进行成像(图17)。可见，无论是成像范围还是成像精度，基于正态分布权函数的VSP-CDP叠加(图17b)的成像结果都明显优于常规VSP-CDP叠加(图17a)。图18为过井地面地震剖面与非零井源距VSP纵波时间域成像剖面。由图可见：非零井源距VSP纵波时间域成像剖面(图18b)与过井地面地震剖面(图18a)波组特征一致，地层产状相近；在2400～2700ms层段后者能量较弱，地层内幕结构信息模糊(图18a)，前者的波组特征连续，成层性好，内幕信息丰富，分辨率高于后者(图18b)。