0 引言

由于地层的吸收衰减作用，地震波在传播过程中会逐渐衰减，地面地震数据通常具有有限带宽，这可能会导致无法成像或无法识别细微的地层特征^[1-3]。影响地面地震数据分辨率的主要因素之一是近地表松散风化带对地震波的强烈吸收，然而，与地面地震勘探不同，VSP是在井中接收地震波，其数据受近地表的影响较小，因此，VSP数据通常具有比地面地震数据更宽的频带以及更高的信噪比和分辨率^[4-6]。目前有4种方法可以通过VSP数据来提高地面地震数据的分辨率：①从VSP数据中提取品质因子Q ^[7-8]，然后将反Q滤波方法应用于叠后地表地震记录^[9-11]；②高频恢复方法^[12-13]，从不同深度的VSP数据中提取滤波算子，然后对叠后地面地震数据进行高频恢复滤波；③VSP小波匹配方法，分别从叠后地面地震数据和VSP数据中提取子波，然后将二者的子波进行匹配^[14-15]，以改善叠后地面地震数据的分辨率；④VSP子波反褶积方法，从VSP数据中提取子波，然后计算反褶积算子，并对叠后地表地震数据进行反褶积^{[14, 16]}。

目前的VSP子波反褶积方法仅从VSP地震记录中提取一个时不变的反褶积算子，但由于地下介质的黏弹性，VSP记录的子波是逐渐衰减的，因为它忽略了子波的时变特性，所以并不准确。为了改进该方法，本次研究将零偏移距VSP数据的下行直达波作为单程地震子波，并将其振幅谱转换为双程地震子波振幅谱，在地震子波是最小相位的假设下，从振幅谱中获得时变反褶积算子，以期其可用于对叠后地面地震数据进行时变子波反褶积。

1 方法原理

本次研究提出的时变子波反褶积方法，首先使用零偏移距VSP资料的下行直达波来获得时变子波振幅谱，以准确反映单程地震子波的振幅谱随旅行时间的变化；然后利用直达波自身的衰减规律将单程地震子波的振幅谱转换为双程地震子波的振幅谱。该方法的步骤为：①利用零偏移距VSP数据估算单程地震子波振幅谱；②基于单程地震子波振幅谱和地层衰减函数，估算叠后地面地震数据的双程时变子波振幅谱；③在最小相位假设下计算时变反褶积算子；④通过地面地震记录的Gabor谱和反褶积算子的乘积得到反射系数的Gabor谱；⑤通过反Gabor变换计算反射系数(图 1)。

1.1 时变子波估算

假设W_vsp(t_n, f)和W_vsp(t_{n + 1}, f)分别为零偏移距VSP记录下行直达波在t_n和t_{n + 1}时刻的振幅谱，且n = 1，2，…，N - 1，N为检波器的个数。根据地层的吸收衰减，它们满足^[17]

$ \ln \left[\frac{W_{\mathrm{vsp}}\left(t_{n+1}, f\right)}{W_{\mathrm{vsp}}\left(t_{n}, f\right)}\right]=\frac{-\pi f\left(t_{n+1}-t_{n}\right)}{Q} $

(1)

式中：t为时间，s；f为频率，Hz；Q为品质因子。

定义c(f) = -π f/Q，Δ t = t_{n + 1} - t_n，则式(1)变为

$ c(f)=\frac{1}{\Delta t\left[\frac{W_{\mathrm{vsp}}\left(t_{n+1}, f\right)}{W_{\mathrm{vsp}}\left(t_{n}, f\right)}\right]} $

(2)

式中：c_n(f)表示t_n和t_{n + 1}间的衰减函数。

对于任意时刻t(t_n < t < t_{n + 1})，VSP单程地震子波的振幅谱W_vsp(t, f)可写为

$ W_{\mathrm{vsp}}(t, f)=W_{\mathrm{vsp}}(t, f) \exp \left[c_{n}(f)\left(t-t_{n}\right)\right] $

(3)

将c₁(f)作为0到t₁时刻的衰减函数，就可以计算0到t_{n + 1}任意时刻的单程地震子波振幅谱。

由于VSP地震子波是单程衰减的振幅谱，因此需要将VSP单程地震子波的振幅谱转化为与地面地震相对应的双程地震子波的振幅谱。对于任意时刻t，当它满足t_n < t < t_{n + 1}时，则有

$ \begin{array}{l} {W_{{\rm{surface }}}}(2t, f) = {W_{{\rm{rsp }}}}(t, f){\rm{ exp}} \cdot \\ {\rm{ }}\left[ {\sum\limits_{m = 1}^{n - 1} {{c_m}} (f){t_m} + {c_n}(f)\left( {t - {t_n}} \right)} \right] \end{array} $

(4)

式中：W_surface(2 t, f)为叠后地面地震双程地震子波的频谱。

1.2 时变地震子波反褶积

信号u(t)的Gabor变换定义为^[18]

$ U(t, f)=\int_{-\infty}^{+\infty} u(\tau) g(\tau-t) e^{-i 2 \pi f t} \mathrm{d} \tau $

(5)

式中：U (t, f)为信号u(τ)的Gabo r谱；g(τ)为Gabor分析时窗；t为分析时窗的中心位置，s。

地震记录的Gabor谱可以近似写为3个部分的乘积^[19-21]，表达式为

$ U(t, f) \approx W(t, f) R(t, f)=W(f) R(t, f) \exp \left(\frac{-\pi f t}{Q}\right) $

(6)

式中：U (t, f)是地震记录u(t)的Gabor谱；W (f)是子波w(t)的傅里叶变换；R(t, f)是反射系数r(t)的Gabor变换。

一旦计算得到时变地震子波谱W_surface(t, f)，在最小相位假设下，反褶积算子可写为^[19]

$ D(t, f) = \left[ {\frac{1}{{{W_{{\rm{surface }}}}(t, f)}}} \right]{e^{ - iH\left[ {{W_{{\rm{surface }}}}(t, f)} \right]}} $

(7)

式中：H[]表示Hilbert变换。

根据式(6)，反射系数R_g(t, f)的Gabor谱可以由下式计算得到

$ R_{\mathrm{g}}(t, f)=U(t, f) D(t, f) $

(8)

再利用Gabor反变换即可计算得到反射系数。

2 理论模型数据试算

图 2(a)为层状模型和VSP观测系统。在该模型中，品质因子Q的数值是逐渐增加的。在100 m到2 080 m深度内布置100个检波器，检波器间距为20 m。震源与井口之间的水平距离为20 m，地震记录长度为1.2 s，采样间隔为0.002 s。由图 2 (b)可以看到，该记录中直达波的振幅和主频随深度减小，表现出单程VSP子波的时变特性。

从图 3可看出，直达波的振幅和主频都随着旅行时的增加而减小。利用式(1)-(3)可计算出每个时间采样点的VSP单程子波振幅谱[图 3 (b)]，利用式(4)，可得到双程地震子波的振幅谱[图 3 (c)]。比较图 3 (b)和图 3 (c)可知，双程地震子波振幅谱的衰减明显强于单程地震子波的振幅谱。为了验证从VSP单程地震子波振幅谱计算的双程地震子波振幅谱的正确性，使用图 2 (a)的层状模型直接计算理论双程地震子波的振幅谱[图 3(d)]。比较图 3(c) 和图 3(d)可以看到，通过计算得到的双程地震子波的振幅谱与理论子波振幅谱几乎一致。

在估算得到双程地震子波振幅谱后，可以对理论合成地震记录进行时变子波反褶积(图 4)。图 4 (b) 中的理论合成衰减地震记录是由图 3 (d)中的时变子波和图 4 (a)中的理论随机反射系数褶积产生。利用式(7)和式(8)来进行时变子波反褶积，其中反褶积算子由图 3(c)中所示的时变子波振幅谱计算得到。由时变子波反褶积后的地震记录[图 4(c)]可以看到，时变子波反褶积后的记录更接近理论反射系数。

3 实际地震数据应用

由图 5 (a)可以看到，下行直达波的振幅能量很强，而上行波的振幅能量较弱，随着传播距离的增加，下行直达波的振幅在逐渐衰减。对零偏移距VSP地震记录进行波场分离，并将下行直达波截断为单程VSP地震子波，地震记录如图 5 (b)所示。

图 6为利用实际VSP数据下行直达波的振幅谱计算得到双程地震子波的振幅谱。由图 6 (a)，可以看到，随着深度和旅行时的增加，振幅减小，频带变窄。图 6(b)和图 6(c)分别为针对每个采样点计算得到的VSP单程子波的幅度谱和双程子波的幅度谱。将二者进行对比可知，双程子波的振幅谱具有比单程子波的振幅谱更弱的振幅和更窄的频带。

在获得双程时变子波振幅谱后，在最小相位假设下计算反褶积算子。由图 7(a)可以看到，随着深度的增加，地震记录的振幅能量减弱，同相轴变粗，分辨率降低。利用上述得到的反褶积算子对其进行时变子波反褶积处理[图 8 (b)]可以看到，时变反褶积后，中深层地震记录的能量得到有效补偿，且同相轴被压缩，波组信息更加丰富，分辨率明显提高。

4 结论

(1) 提出一种基于零偏移距VSP的时变子波反褶积方法，该方法利用零偏移距VSP数据的直达波估算双程时变子波的振幅谱，然后利用双程时变子波的振幅谱估算得到时变子波的反褶积算子，将其用于叠后地面地震记录，其中时变地震子波的估算精度会影响最终的反褶积效果。

(2) 理论合成数据和实际地震资料的处理结果证明本文方法可以有效恢复地震记录的高频衰减能量，并可提高地震记录的分辨率。