地层倾角信息对于地震采集、处理、解释都具有重要意义。在采集阶段, 地层倾角信息可以用于指导观测系统设计; 在处理阶段, 地层倾角对于速度分析、速度建模等有很大影响; 在油藏描述阶段, 真实的地层倾角信息有助于研究地质构造和沉积环境, 进而追踪地下油气的分布情况。地层倾角的计算方法有局部倾斜叠加法^[1]、复地震道分析法^[2]、结构张量法^[3-4]、平面波分解法^[5]以及相干估计法^[6-7]等。CLAERBOUT^[5]提出了平面波分解滤波器, 通过有限差分法求解平面波微分方程得到地震数据的局部斜率。一般情况下, 局部平面波模型能够方便地表达地震数据, 但是需要预测较多的滤波器系数, 而这些系数需要人工交互处理, 存在一定的局限性。FOMEL^[8]对CLAERBOUT滤波器^[5]进行了改进, 结合适当的数据正则化方法, 使滤波器系数预测减少到最少, 局部平面波斜率预测结果具有清晰的物理意义。FOMEL^[9]将改进的平面波分解滤波器应用于局部斜率估计、断层预测、数据插值和噪声衰减。杜婧等^[10]使用平面波分解滤波器估计局部斜率属性, 用于VSP波场分离研究。刘洋等^[11]使用平面波分解滤波器得到地震数据的局部倾角属性, 设计出预测构造数据体和构造走向的滤波器, 该滤波器能在去噪的同时保护断层信息。刘海燕^[12]利用平面波分解滤波器估计的地层倾角对地震数据进行调和, 使其C3相干体算法的分辨率更高, 有效消除了地层倾角导致的低相干干扰, 突出了有效断层信息。陈鑫^[13]在速度无关时差校正中使用平面波分解滤波器获得局部同相轴斜率信息, 用于替代相关速度参数进行动校正处理。

分析发现, 常规平面波分解方法只利用地震数据, 未充分利用老油区的大量倾角测井数据。本文基于平面波分解滤波器计算局部斜率的方法, 将地层倾角测井信息作为先验信息, 结合Tikhonov正则化思想, 提出了井震融合平面波分解地层倾角计算方法; 利用理论数据和实际资料对方法进行了测试。

1 平面波分解方法原理

根据局部平面波的物理模型, 可以定义局部平面波微分方程:

$ \frac{{\partial P\left( {t,x} \right)}}{{\partial x}} + \sigma \frac{{\partial P\left( {t,x} \right)}}{{\partial t}} = 0 $

(1)

式中:P(t, x)代表局部平面波波场, σ代表局部斜率, 两者都与时间t和炮检位置x有关。公式(1)表明, 在x₁处观测到的波场可以被x₂处观测到的波场完美地预测, 它们相差一个延迟时间σΔx, 其中Δx=|x₁-x₂|。即平面波时延满足:

$ P\left( {t,x} \right) = P\left( {t + \sigma \Delta x,x + \Delta x} \right) $

(2)

将方程(2)对时间和空间分别进行Z变换, 得:

$ \left( {1 - Z_t^\sigma {Z_x}} \right)P\left( {{Z_t},{Z_x}} \right) = 0 $

(3)

式中:P(Z_t, Z_x)为P(t, x)的Z变换, 1-Z_t^σZ_x为平面波分解算子。利用一个全通滤波器$\frac{{B({Z_t})}}{{B(1/{Z_t})}}$近似时移算子Z_t^σ, 平面波分解算子1-Z_t^σZ_x可以表示为^[7]:

$ C\left( {Z_t^\sigma ,{Z_x}} \right) = B\left( {\frac{1}{{{Z_t}}}} \right) - {Z_x}B\left( {{Z_t}} \right) $

(4)

其中,

$ B\left( {{Z_t}} \right) = \sum\limits_{k = - N}^N {{b_k}\left( \sigma \right)Z_t^{ - k}} $

(5)

式中:N为空间滤波器的阶数; b_k(σ)为局部斜率σ的函数。在低频处(N=1), 通常采用三点中心滤波器, 即:

$ \begin{array}{l} {B_3}\left( {{Z_t}} \right) = \frac{{\left( {1 - \sigma } \right)\left( {2 - \sigma } \right)}}{{12}}Z_t^{ - 1} + \frac{{\left( {2 + \sigma } \right)\left( {2 - \sigma } \right)}}{6} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{\left( {1 + \sigma } \right)\left( {2 + \sigma } \right)}}{{12}}{Z_t} \end{array} $

(6)

方程(4)相当于一个二维滤波算子, 将方程(4)代入方程(3), 可以得到与地层倾角σ相关的方程:

$ C\left( {Z_t^\sigma ,{Z_x}} \right)P\left( {{Z_t},{Z_x}} \right) \approx 0 $

(7)

方程(7)即利用平面波分解滤波器计算地层局部斜率的基本公式。

2 井震融合的地层倾角计算

在方程(7)中, 我们重点关注局部斜率σ。为了与最优化理论常用公式一致, 将数据P(Z_t, Z_x)简化为d, C(Z_t^σ, Z_x)简化为C(σ), 得:

$ C\left( \sigma \right)d \approx 0 $

(8)

方程(8)的本质是求解一个非线性问题。根据最优化理论, 可将方程(8)转化为线性方程(9), 然后使用高斯-牛顿法求解:

$ C'\left( \sigma \right)\Delta \sigma d + C\left( \sigma \right)d \approx 0 $

(9)

式中:Δσ表示斜率的增量, C′(σ)是滤波器C(σ)的导数。在迭代求解线性方程(9)的过程中, 初始σ值被不断更新, 逐步趋近于真实斜率。

在地震数据信噪比较低、分辨率不高的情况下, 直接求解方程(9)往往不稳定, 需要引入合适的正则化。FOMEL^[7]给出正则化项:

$ \varepsilon D\Delta \sigma \approx 0 $

(10)

式中:D为梯度算子, ε为正则化比例系数。通过增加方程(10)中正则化项, 可以得到更加可靠的局部斜率。

在成熟探区, 通常有各种测井数据, 由井口附近的倾角测井数据能够得到较地震数据更为可靠的局部倾角信息。根据Tikhonov正则化的思想, 在公式(8)中加入先验测井信息约束, 则得到井震融合的倾角估计方程:

$ C\left( {\sigma + \xi \sigma } \right)d \approx 0 $

(11)

式中:ξσ表示引入的测井约束, 本文主要利用地层倾角测井信息。

由公式(11)得到的是沿主测线和联络测线方向的局部倾角, 称之为视倾角, 需要经过一定的三角换算转化为真倾角和真倾向。转换公式如下:

$ \tan \theta \cdot \cos \varphi = \tan {\theta _y} $

(12a)

$ \tan \theta \cdot \sin \varphi = \tan {\theta _x} $

(12b)

$ \varphi = \arctan \left( {\frac{{\tan {\theta _x}}}{{\tan {\theta _y}}}} \right) $

(13a)

$ \theta = \arctan \left\{ {\frac{{\tan {\theta _x}}}{{\sin \left\{ {\arctan \left( {\frac{{\tan {\theta _x}}}{{\tan {\theta _y}}}} \right)} \right\}}}} \right\} $

(13b)

式中:θ_x和θ_y分别为沿主测线和联络测线的视倾角; θ为地层真倾角; φ为真倾向与y方向视倾向之间的夹角(正北方为y方向, 正东方为x方向)。由公式(13)可以得到真实的地层倾角和倾向。

基于地层倾角信息约束的平面波分解地层倾角估计方法有两方面优势:①将经过预处理后的倾角测井信息作为约束加入到方程(11)中, 有倾角测井值的网格不更新σ, 在一定程度上减少了反演解的误差, 降低了反演迭代次数, 提高了反演效率。②利用倾角测井信息约束提高地层倾角估计的精度, 有助于精细刻画薄互层、微幅构造等小尺度地质体。

3 测试与应用 3.1 理论模型测试

采用图 1a所示合成3D偏移数据体^[6]对本文方法进行了测试, 图 1b为合成数据的平面图以及4口合成地层倾角测井(垂直井)井口位置示意图。将图 1b中红五星周围3口井作为约束井, 采用公式(11)得到x方向和y方向的视倾角如图 2所示。再利用公式(13)求出地层真倾角及倾向, 如图 3所示。为了进一步说明本文方法的正确性, 我们对比了平面波分解方法未加入约束、加入井约束预测的倾角曲线以及合成地层倾角测井曲线, 如图 4所示。可以看出, 加入井约束预测的倾角更加趋近于地层倾角测井的结果。

3.2 实际资料应用 3.2.1 地层倾角估计

南堡凹陷为经过多期构造运动形成的裂陷盆地, 地质特征复杂。位于南堡凹陷南部的南堡3号(NP3)构造地层倾角变化大, 严重影响地质认识程度。采用本文井震融合的地层倾角估计方法估算了南堡3号构造的地层倾角数据体。图 5a为NP3号构造1599线偏移剖面, 利用公式(11)迭代计算了地层视倾角/向, 图 5b和图 5c为利用公式(13)由视倾角/向转换的地层真倾角/向。

根据南堡3号工区NP3-19井、NP306X1井和PG2井的轨迹, 抽取未加入井约束地层倾角估计方法得到的地层倾角数据体中相应的倾角值, 将其与地层倾角测井结果作对比, 如图 6所示。

由于地层倾角测井的采样率和井轨迹数据采样率不一致, 故两者不可能一一对应, 但是从图 6可以清晰地看出, 未加入井约束地层倾角估计方法得到的3口井地层倾角变化趋势与地层倾角测井结果基本吻合。

将NP3-19井、NP306X1井和PG2井测井数据作为约束加入公式(11)进行地层倾角估算, 再抽取离NP3-19井较近的NP3-20井地层倾角信息进行对比, 结果如图 7所示。由图 7可见, 加入测井数据约束后估算出来的地层倾角更加趋近于地层倾角测井的结果, 估算精度得到了提高。

对比PG2, NP3-20以及P3002三口井的地层产状(表 1)可以看出, 加入测井数据约束后, 地震数据体地层产状估算精度明显提升, 与实际测量结果相对误差小于1%。P3002井与约束井距离较远(约800m左右), 约束效果不明显, 地层产状相对误差为1%左右, 也在可接受范围内。

3.2.2 地层厚度校正

为了说明井震融合平面波分解地层倾角估计在薄互层、微幅构造等小尺度地质体精细刻画中的作用, 这里以地层厚度校正为例, 对研究区内NP3-82井、PG2井、NP3-80井进行测井地层厚度校正, 将地层视厚度换算成地层真厚度。

根据图 5所示的地层倾角信息, 利用测井曲线校正公式^[14], 将NP3-82井、PG2井、NP3-80井的垂直厚度转换成真厚度(图 8)。图 8a为NP3-82井-PG2井-NP3-80井地层垂直厚度对比图; 图 8b和图 8c分别为未加入井约束和加入井约束方法估算地层倾角后, 校正得到的NP3-82井-PG2井-NP3-80井地层真厚度对比图。从图 8可以看出, PG2井和NP3-80井地层水平厚度变化较大, 与海相地层沉积稳定的地质规律不符(图 8a); 采用未加入井约束方法估算地层倾角并校正地层厚度后, PG2井和NP3-80井地层水平厚度变化趋势有所减缓(图 8b); 而采用本文井震融合方法预测地层倾角并对地层厚度进行校正后, 地层水平厚度变化趋势更加平稳, 更加符合海相地层地质规律(图 8c), 小层划分精度也得到了提高。

4 结束语

本文提出一种井震融合的平面波分解地层倾角估算方法, 根据Tikhonov正则化思想, 融合地震和地层倾角测井信息实现地层倾角的估计。该方法引入测井信息约束后, 不仅使平面波分解滤波器估计地层倾角的精度得到了提升, 减小了与实测结果的误差, 而且提高了迭代反演的收敛速率, 因此具有一定的实际应用价值。