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 动态射线追踪

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

$ \frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial z^{2}}=\frac{1}{v} \frac{\partial^{2} u}{\partial t^{2}} $

(1)

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

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

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

$ \frac{1}{h} u_{, ss}+h u_{, nn}-\frac{h}{v^{2}} u_{, tt}+h u_{, s}\left(\frac{1}{h}\right)_{, s}+u_{, n} h_{, n}=0 $

(2)

式中：$ h=1+\frac{v, n}{v} n$是坐标变换中的标度因子，其中v=v(s, n)为地震波速度，$v_{, n}=\left[\frac{\partial v(s, n)}{\partial n}\right]_{n=0} ; u_{, s s} $、u_{, nn}、u_{, tt}分别为u对s、n、t的二阶偏导数；u_{, s}、u_{, n}分别为u对s、n的一阶偏导数；$ \left(\frac{1}{h}\right)_{s}$为1h对s的一阶偏导数；h_{, n}为h对n的一阶偏导数。

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

$ u(s, n, t)=\exp \left\{-\mathrm{i} \omega\left[t-\int_{s_{0}}^{s} \frac{\mathrm{d} s}{v(s)}\right]\right\} U(s, n, \omega) $

(3)

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

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

$ u(s, n, t)=\sqrt{v(s)} \exp \left\{-\mathrm{i} \omega\left[t-\int_{s_{0}}^{s} \frac{\mathrm{d} s}{v(s)}\right]\right\} W(s, v) $

(4)

其中

$ W(s, v)=A(s) \exp \left(\frac{\mathrm{i}}{2} v^{2} \varGamma\right) $

(5)

$ \varGamma_{, s}+v \varGamma^{2}+v^{-2} v_{, nn}=0 $

(6)

$ A_{, s}+\frac{1}{2} v A \varGamma=0 $

(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)，设

$ \varGamma(s)=\frac{1}{v q} q_{, s} $

(8)

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

$ v q_{, ss}-v_{, s} q_{, s}+v_{, nn} q=0 $

(9)

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

$ v(v_{, s} q+v p_{, s})-v_{, s}(v p)+v_{, nn} q=0 $

(10)

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

$ \left\{\begin{array}{l} q_{, s}=v p \\ p_{, s}=-v^{-2} v_{, m} q \end{array}\right. $

(11)

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

$ \varGamma=p q^{-1} $

(12)

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

$ \begin{aligned} u(s, n, t)=& \varPsi\left[\frac{v(s)}{q(s)}\right]^{\frac{1}{2}} \exp \left\{-\mathrm{i} \omega\left[t-\int_{s_{0}}^{s} \frac{\mathrm{d} s}{v(s)}\right]+\right.\\ &\left.\frac{\mathrm{i} \omega}{2} \frac{p(s)}{q(s)} n^{2}\right\} \end{aligned} $

(13)

式中$ \varPsi\left[\frac{v(s)}{q(s)}\right]$为权函数。引入如下记号

$ \left\{\begin{array}{l} K(s)=v(s) \operatorname{Re}\left[\frac{p(s)}{q(s)}\right] \\ L_{\mathrm{h}}=\left\{\frac{\omega}{2} \operatorname{Im}\left[\frac{p(s)}{q(s)}\right]\right\}^{-\frac{1}{2}} \\ \tau(s)=\int_{s_{0}}^{s} \frac{\mathrm{d} s}{v(s)} \end{array}\right. $

(14)

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

$ \begin{aligned} u(s, n, t)=& \varPsi\left[\frac{v(s)}{q(s)}\right]^{\frac{1}{2}} \exp \{-\mathrm{i} \omega[t-\tau(s)]+\\ &\left.\mathrm{i} \omega \frac{n^{2}}{2 v(s)} K(s)-\frac{n^{2}}{L_{\mathrm{h}}^{2}}\right\} \end{aligned} $

(15)

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

$ U(s, t)=\int_{0}^{2 \pi} \varPsi(\varphi) u_{\varphi}(s, n) \mathrm{d} \varphi $

(16)

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

1.2 正态分布权函数

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

$ A=\cos ^{2}\left(\frac{\pi}{2} n_{s}^{2} L_{\mathrm{h}}^{-2}\right) $

(17)

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

$ f(x)=\frac{1}{\sqrt{2 \pi} \sigma} \mathrm{e}^{-\frac{(x-\mu)^{2}}{2 \sigma^{2}}} \quad-\infty<x<+\infty $

(18)

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

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

$ A=\frac{1}{\sqrt{2 \pi} L_{\mathrm{h}}} \mathrm{e}^{-\frac{n_{{s}}^{2}}{2 L_{\mathrm{h}}^{2}}} $

(19)

当相邻射线与中心射线重合时(n_s=0)，由单一中心射线的波场能量表征该接收道的波场(常规射线追踪方法)；当n_s < L_h时，由L_h范围内所有射线波场能量的加权表征该接收道的波场(旁轴射线追踪方法)；当n_s>L_h时，相邻射线的波场能量对该接收道的波场无贡献。

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

1.3 VSP-CDP转换 1.3.1 常规VSP-CDP转换原理

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

$ t=\frac{\sqrt{x_{0}^{2}+\left(2 h-z_{\mathrm{R}}\right)^{2}}}{V} $

(20)

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

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

$ \left\{\begin{array}{l} x_{1}=\frac{h x_{0}}{2 h-z_{\mathrm{R}}} \\ z_{1}=h \end{array}\right. $

(21)

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

图 5为单道数据的VSP-CDP转换示意图。由图可见，根据上述VSP-CDP转换方法可以计算到

达检波点R的两个反射界面的反射点坐标(x₁, z₁)和(x₂, z₂)，分别对应T₁和T₂时刻，VSP-CDP转换将共炮点道集记录中T₁、T₂时刻之间的地震记录分别归位到反射点(x₁, z₁)、(x₂, z₂)处(若要转换成井源距—反射点深度域，需要对T₁时刻之前的地震记录进行时深转换，再进行反射点归位)，这样就形成了两个反射界面情况下VSP单道数据的VSP-CDP转换。同理，可完成单道数据所有反射界面及整个VSP多道数据的VSP-CDP转换。

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

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

2 模型试算

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

3 方法对比

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

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)。

5 结束语

基于动态射线追踪有效邻域波场近似理论，本文提出了一种非零井源距VSP资料叠加成像方法。通过研究正态分布的性质及特点，导出了基于正态分布叠加的权函数计算公式，并用于VSP叠加成像，改善了反射点覆盖次数不均匀现象，提高了地震成像横向分辨率。

模型试算及实际资料处理结果表明，基于正态分布权函数的VSP-CDP叠加成像方法能够对复杂构造VSP资料精确成像，验证了方法的有效性和稳定性。