0 引言

近年来，准噶尔、松辽和渤海湾等多个含油气盆地的火山岩油气藏已经成为重要的油气增储上产领域^[1]。由于火山机构的改造强烈、蚀变严重、裂缝发育、岩性多变和多期叠置，导致火山岩的地震响应特征不明确。准确识别火山机构是火山岩气藏勘探面临的一个突出问题，直接影响对火山岩岩相和火山岩优质储层空间分布的认识，进而影响火山岩油气勘探成效。多年来，国际上更偏重于研究火山岩地层地质特征及成因条件。20世纪70年代末，Humter等^[2]将火山分为单一型、复合型、盾形和古潜火山等。Pike^[3]根据火山机构的形状、成分和成因等对火山分类，按照发育规模，将层状火山分为小型、中型和大型，将盾形火山分为小型和大型。依据火山岩成因条件和火山机构外貌形态等因素，科普切弗—德沃尔尼科夫等^[4]将火山岩岩相分为喷发相、通道相及次火山岩相等。Fisher等^[5]研究了火山碎屑岩的沉积条件及搬运方式，并将其岩相分为火山灰流、碎屑流和喷发冲击相等。在中国，由于完全不同的地质成因，火山岩在很多方面与沉积岩差别巨大。作为沉积事件的产物，火山岩内单元形态多样、规模差异很大、横向迁移非均质性很强、岩石类型多、岩性岩相复杂，致使火山岩具有独特的地震反射特征。基于火山岩储层地质学理论，王璞珺等^[6-7]根据徐家围子断陷岩石露头和钻井信息，将火山岩分为四种相及其十三种亚相，并指出了相识别标志。潘建国等^[8]研究了准噶尔盆地古生界火山岩油气藏的地质特征、物理特性以及成藏规律等。唐华风等^[9]利用自组织神经网络方法对波形分类，采用多参数多次迭代训练、划分模型道，结合单井相分类结果划分松辽盆地营城组上部火山岩平面地震相，结果与地质背景相吻合。

对于火山机构的刻画，大多还是以地震属性方法为主，且都基于全叠加数据^[10-14]。随着地质勘探目标趋于复杂，对地震勘探精度的要求也越来越高，宽方位地震勘探已成为现阶段的主要发展方向之一。对于油气勘探，无论是构造勘探、岩性勘探还是油藏流体检测，宽方位地震技术都是必备的手段，也是在不同阶段刻画地质目标的基础^[15]。五维地震数据可以更好地分析地震波在各向异性介质中传播的旅行时间、速度、振幅、频率和相位等属性随方位角的变化，而且炮检距信息与目标地质体的尺度、地层岩性和流体成分等存在相关性，方位角信息则与地层中的断裂、裂缝等的发育特征相关。

为了降低火山岩预测的多解性，本文利用五维地震资料刻画火山机构，为高效勘探、开发火山岩储层奠定了基础。

1 五维地震数据刻画火山机构

火山机构是指一定时间范围内，由同喷发源的火山物质围绕源区堆积而成的具有一定形态和共生组合关系的各种火山作用产物的总和。目前，大多利用地震剖面的振幅包络面与地层不整合面刻画火山机构的宏观展布特征。由于火山口顶部为喷发能量异常高值区，火山通道内呈高频、低连续、杂乱反射特征，火山通道外地震波速度及能量相对于围岩有较大突变，导致叠后地震属性突变，因此可以利用地震多属性识别火山机构^[16-17]。业界主要基于相干体属性预测火山机构，即利用相邻地震信号的相似性描述地层和岩性的横向不均匀性。另外，由于火山岩埋藏较深，导致基质多变、岩性复杂、岩相变化大。因此，整体上深层地震资料信噪比较低，火山岩地震响应特征不明确、火山机构刻画及火山岩储层预测难度大。当地震波遇到地下物性不连续面时会引起方位各向异性，利用垂直断层走向的地震数据识别断裂构造的效果最佳^[18]，因此本文提出基于五维地震资料的多方位相干算法。

1.1 倾角成像增强处理

在地震属性分析之前，首先对地震数据预处理，包括：对原始地震资料的多尺度、多方位倾角成像增强；利用倾角及方位角的变化计算相邻道的相似性，以提高地震资料横向信噪比。首先，定义三维地震数据体的分方位、分尺度特征$ {\psi }_{u,v}\left(\boldsymbol{x}\right) $，即

$ {\psi }_{u,v}\left(\boldsymbol{x}\right)=\frac{\left|\right|{\boldsymbol{k}}_{u.v}\left|\right|}{{\sigma }^{2}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\left|\right|{\boldsymbol{k}}_{u,v}\left|{|}^{2}\right|\left|\boldsymbol{x}\right|{|}^{2}}{2{\sigma }^{2}}\right)\times \\\;\;\;\;\left[\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{i}{\boldsymbol{k}}_{u,v}\boldsymbol{x}\right)-\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\sigma }^{2}}{2}\right)\right] $

(1)

式中：$ \boldsymbol{x}=(a,b,c)\mathrm{为}\mathrm{空}\mathrm{间}\mathrm{任}\mathrm{意}\mathrm{点},a\mathrm{、}b\mathrm{、}c $$ \mathrm{分}\mathrm{别}\mathrm{代}\mathrm{表} $$ \mathrm{I}\mathrm{n}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{、}\mathrm{C}\mathrm{r}\mathrm{o}\mathrm{s}\mathrm{s}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{、}\mathrm{时}\mathrm{间}\mathrm{方}\mathrm{向}\mathrm{的}\mathrm{尺}\mathrm{度};\sigma $为滤波器带宽；$ {\boldsymbol{k}}_{u,v}={k}_{v}{\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{s}\mathrm{i}\mathrm{n}\theta \right)}^{\mathrm{T}} $，$ {k}_{v}={2}^{\left(2-v\right)/2}\mathrm{\pi } $，$ \theta =u\mathrm{\pi }/K $，u为滤波器的方向参数，决定$ \theta $角的大小，K为核函数，v为滤波器的尺度。

对于给定的任意一个空间点x的三维地震数据$ I\left(\boldsymbol{x}\right), $定义不同方位、不同尺度的地震数据$ {D}_{u,v}\left(\boldsymbol{x}\right) $

$ {D}_{u,v}\left(\boldsymbol{x}\right)={\psi }_{u,v}\left(\boldsymbol{x}\right)*I\left(\boldsymbol{x}\right) $

(2)

对D进行重新编码，得到$ {D}_{u,v}\left(\boldsymbol{x}\right) $在多尺度、多方位滤波数据的二进制编码，并计算合成数据，其实部和虚部分别为$ \boldsymbol{G} $^Re和$ \boldsymbol{G} $^Im。根据分解数据定义张量矩阵，将$ {\boldsymbol{G}}_{}^{\mathrm{R}\mathrm{e}}\mathrm{和} $$ \boldsymbol{G} $^Im用特征值和特征向量表示为

$ {\boldsymbol{G}}_{}^{\mathrm{R}\mathrm{e}}=\left({\boldsymbol{S}}_{1}^{\mathrm{R}\mathrm{e}}{\boldsymbol{S}}_{2}^{\mathrm{R}\mathrm{e}}{\boldsymbol{S}}_{3}^{\mathrm{R}\mathrm{e}}\right)\left(\begin{array}{ccc}{\lambda }_{1}^{\mathrm{R}\mathrm{e}}& 0& 0\\ 0& {\lambda }_{2}^{\mathrm{R}\mathrm{e}}& 0\\ 0& 0& {\lambda }_{3}^{\mathrm{R}\mathrm{e}}\end{array}\right){\left({\boldsymbol{S}}_{1}^{\mathrm{R}\mathrm{e}}{\boldsymbol{S}}_{2}^{\mathrm{R}\mathrm{e}}{\boldsymbol{S}}_{3}^{\mathrm{R}\mathrm{e}}\right)}^{\mathrm{T}} $

(3)

$ {\boldsymbol{G}}_{}^{\mathrm{I}\mathrm{m}}=\left({\boldsymbol{S}}_{1}^{\mathrm{I}\mathrm{m}}{\boldsymbol{S}}_{2}^{\mathrm{I}\mathrm{m}}{\boldsymbol{S}}_{3}^{\mathrm{I}\mathrm{m}}\right)\left(\begin{array}{ccc}{\lambda }_{1}^{\mathrm{I}\mathrm{m}}& 0& 0\\ 0& {\lambda }_{2}^{\mathrm{I}\mathrm{m}}& 0\\ 0& 0& {\lambda }_{3}^{\mathrm{I}\mathrm{m}}\end{array}\right){\left({\boldsymbol{S}}_{1}^{\mathrm{I}\mathrm{m}}{\boldsymbol{S}}_{2}^{\mathrm{I}\mathrm{m}}{\boldsymbol{S}}_{3}^{\mathrm{I}\mathrm{m}}\right)}^{\mathrm{T}} $

(4)

式中：$ {\lambda }_{1}^{\mathrm{R}\mathrm{e}} $、$ {\lambda }_{2}^{\mathrm{R}\mathrm{e}} $、$ {\lambda }_{3}^{\mathrm{R}\mathrm{e}} $为$ {\boldsymbol{G}}_{}^{\mathrm{R}\mathrm{e}} $的三个非负特征值；$ {\boldsymbol{S}}_{1}^{\mathrm{R}\mathrm{e}} $、$ {\boldsymbol{S}}_{2}^{\mathrm{R}\mathrm{e}} $、$ {\boldsymbol{S}}_{3}^{\mathrm{R}\mathrm{e}} $分别为$ {\lambda }_{1}^{\mathrm{R}\mathrm{e}} $、$ {\lambda }_{2}^{\mathrm{R}\mathrm{e}} $、$ {\lambda }_{3}^{\mathrm{R}\mathrm{e}} $对应的特征向量；$ {\lambda }_{1}^{\mathrm{I}\mathrm{m}} $、$ {\lambda }_{2}^{\mathrm{I}\mathrm{m}} $、$ {\lambda }_{3}^{\mathrm{I}\mathrm{m}} $为$ {\boldsymbol{G}}_{}^{\mathrm{I}\mathrm{m}} $的三个非负特征值，且满足$ {\lambda }_{1}^{\mathrm{I}\mathrm{m}} $≥$ {\lambda }_{2}^{\mathrm{I}\mathrm{m}} $≥$ {\lambda }_{3}^{\mathrm{I}\mathrm{m}} $≥$ 0 $；$ {\boldsymbol{S}}_{1}^{\mathrm{I}\mathrm{m}} $、$ {\boldsymbol{S}}_{2}^{\mathrm{I}\mathrm{m}} $、$ {\boldsymbol{S}}_{3}^{\mathrm{I}\mathrm{m}} $分别为$ {\lambda }_{1}^{\mathrm{I}\mathrm{m}} $、$ {\lambda }_{2}^{\mathrm{I}\mathrm{m}} $、$ {\lambda }_{3}^{\mathrm{I}\mathrm{m}} $对应的特征向量。

基于特征值和特征向量求取倾角属性，利用倾角及方位角变化计算相邻道的相似性，以提高地震资料垂向分辨率，并保留其不连续的波组特征，使地震数据平滑，从而更容易识别不连续性(如断层)。图 1为原始地震数据及其倾角成像增强处理剖面。由图可见，在分方位角倾角成像增强处理剖面上火山机构及其宏观展布特征更清晰，火山机构边界也更清楚（图 1b、图 1d、图 1f）。

1.2 方位时窗分析

通过对不同方位的地震数据，沿与观测方位一致的方向计算相干值，以突出垂直于观测方位的断裂特征。在实际分方位地震数据处理中，通常按30°或22.5°采样间隔划分方位角，一般分为6~8个方位数据体，而常规尺寸的三维分析时窗只能取4个方向。要实现多方位相干，需扩大窗口的平面尺寸，从而提高计算稳定性。通过设置不同分析窗的平面尺寸，增强相干不连续性，使断裂特征更明显。为避免大时窗分析引起分辨率下降，从而使部分信息模糊，因此在方位时窗分析中增加反距离加权，即在计算相干属性时对不同地震道赋予不同权重——增加距中心道近的地震道的权重（减小距中心道远的地震道的权重），从而有效降低分析窗尺寸增大带来的平均效应。通过在同一个方向计算地震道的加权相干，作为该方向中心道的相干值。

设某个方位分析窗包含Y（Y > 5）道（图 2），除中心道及相邻两道外，其余道使用高斯核函数w_y

$ {w}_{y}=\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\left(y-\frac{Y+1}{2}\right)}^{2}}{2{\varphi }^{2}}\right] $

(5)

作为权值函数分别计算每一道的权重，并通过

$ {w}_{y}^{\text{'}}=\frac{2}{\sum \limits_{y=1}^{y}{w}_{y}}{w}_{y}y\ne \frac{Y-1}{2},y\ne \frac{Y+1}{2} $

(6)

确保每一侧的权重之和为1。式中：$ {w}_{y}^{\text{'}} $为第y（y=1, 2, …, Y）道的权重；$ \varphi $为权重因子，$ \varphi $越小表示距中心道距离增加、权重减小，反之亦然。

对中心道两侧的振幅分别加权求和

$ \left\{\begin{array}{l}{d}_{1}=\sum\limits _{y=1}^{\frac{Y-3}{2}}{s}_{y}{w}_{y}^{\text{'}}\\ {d}_{2}=\sum \limits_{y=\frac{Y+5}{2}}^{Y}{s}_{y}{w}_{y}^{\text{'}}\end{array}\right. $

(7)

式中：$ {d}_{1} $、$ {d}_{2} $表示中心道两侧地震道的振幅加权平均；$ {s}_{y} $和$ {w}_{y}^{\text{'}} $分别为地震道和相应的权重值。在大时窗分析时，将$ {d}_{1} $、$ {d}_{2} $作为两端新地震道的振幅，考虑到不同位置地震道的贡献，保持方位分析窗内为5道，可以获得4个不同方向的数据。

图 3为插值后方位分析窗示意图。通过对地震道进行反距离加权插值，在分析窗内可以获得8个不同的方向，并且每个方位分析窗的地震道数相同。

1.3 组合属性技术

由于火山喷发过程中，火山通道及附近会出现塌陷或者崩塌现象，表现为杂乱地震反射。利用五维地震资料，基于倾角成像增强处理计算相干属性再计算梯度模属性，利用分方位属性组合刻画火山机构可反映火山口的细节变化。不同方位梯度模属性火山口预测结果（图 4）表明，火山口主要沿断裂带发育，呈杂乱反射、环形结构，与实钻结果吻合。组合属性技术大大提高了火山机构预测精度，与常规预测方法相比，得到的隐蔽火山口的形态更清晰、火山通道位置更准确，为火山岩岩相预测提供了可靠的数据支撑。

2 属性融合岩相综合预测

在宏观火山机构预测基础上，进一步预测火山岩优势岩相。为预测研究区火山岩岩相，首先利用测井与地震的相互关系建立火山岩喷发模式；在明确不同火山岩岩相的地震响应特征后，再利用方位属性融合法落实火山岩岩相带的平面分布，同时利用“相变法”在地震剖面上精细刻画异常火山岩岩体；最终落实石炭系有利岩体目标。本文主要采用核主成分分析（Kernel Principal Component Analysis，KPCA）融合方位优势属性^[19]，即将原始数据通过非线性映射函数映射到高维特征空间，在高维空间进行主成分分析（Principal Component Analysis，PCA）。使用KPCA融合多属性时，根据核函数首先计算核矩阵K，然后特征分解K，获取最大特征值对应的特征向量α，进而得到权重向量P

$ \boldsymbol{P}= \boldsymbol{\phi } (\boldsymbol{C})\boldsymbol{\alpha} $

(8)

式中$ \boldsymbol{\phi }\left(\boldsymbol{C}\right) $为非线性映射函数，C为由同一道、不同方位的相干值组成的矩阵。

取$ \boldsymbol{\alpha } $与方位属性进行运算，便得到包含大部分信息的方位融合属性。方位相干属性融合过程可表示为

$ {\boldsymbol{C}}_{\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{t}\mathrm{i}}={\boldsymbol{\phi }}^{\mathrm{T}}\left(\boldsymbol{C}\right)\boldsymbol{\phi }\left(\boldsymbol{C}\right)\boldsymbol{\alpha }=\boldsymbol{K}\boldsymbol{\alpha } $

(9)

式中C_multi即为本文定义的多方位相干矩阵。

PCA（图 5a）和KPCA（图 5b）均能识别火山机构，但KPCA方法突出刻画了局部断裂，更利于识别火山机构岩相的平面特征。

综上所述，基于五维地震数据的火山机构刻画及岩相综合识别的技术流程为：首先计算方位几何属性；考虑地下地质体在垂直观测系统方向上响应更明显的优势，构建方位分析时窗提取优势方位信息，利用倾角增强处理地震资料刻画火山机构，明确火山机构宏观分布范围；再结合KPCA融合优势属性预测火山岩有利岩相，精细预测火山岩有利区带（图 6）。

3 应用实例

准噶尔盆地DX地区火山岩具有纵向多期喷发、横向有利火山岩岩相叠置连片的特征。其中，最有利的爆发相火山岩储层主要以“丘状”或“块状”低频、中—弱振幅、弱连续反射为主，溢流相火山岩储层以中—低频、强振幅连续反射为主。火山沉积相则是经常与火山岩共生的一种岩相，可出现在火山活动的各个时期，碎屑成分中含有大量的火山岩岩屑，主要为火山岩穹隆之间的碎屑沉积体，具韵律层理、水平层理。采用单一振幅属性刻画火山岩岩相通常造成不同岩相边界不清楚；采用优势属性融合综合预测火山岩岩相，得到的岩相宏观特征、边界清晰（图 7），利用多属性融合结果可以划分火山岩层序，从而为火山岩储层预测奠定基础。利用多属性融合方法得到了准噶尔盆地KM1井区火山岩岩相分布结果，最终了解了相带宏观展布规律，与井数据符合率达92%（图 8）。

4 结束语

五维地震资料可以更好地分析地震波在各向异性介质中传播的旅行时间、速度、振幅、频率和相位等属性随方位角的变化，本文提出了基于五维地震资料的火山机构刻画和岩相识别技术，获得以下认识：

（1）对原始地震资料进行多尺度、多方位倾角成像增强处理，利用倾角及方位角的变化计算相邻道的相似性，可以提高地震资料的横向信噪比，为火山机构刻画提供了有利数据支撑；

（2）通过定义方位时窗，结合地震道反距离加权算法，提取每个方位对断裂最敏感的信息，提高了火山机构刻画精度，得到的火山形态更清晰；

（3）利用多属性融合预测火山岩岩相避免了单一属性预测火山岩岩性精度低的问题，且多属性融合结果能够有效指导火山岩层序划分，为火山岩储层预测奠定了基础。