0 引言

在地震层序划分的基础上，地震相分析利用地震波参数(振幅、频率、相位等)特征的差异性，将地震层序划分为不同的地震相类别。结合钻井及已知地质信息，可以将地震相分析结果用于岩相或沉积环境的推断^[1]。特定情况下，地震相分析结果也可用于识别储层和流体性质。在地震储层特征的描述和检测方法中，地震相分析已经成为一种不可或缺的方法^[2]。诸多地球物理学学家借鉴图像模式识别方法理论，仅仅从海量地震数据中得到地震属性以及其他一些辅助信息刻画地质体的空间分布，逐渐形成了方法理论较完善的无监督地震相分析方法，实际资料应用也取得了良好效果^[3]。

自组织映射(Self-Organizing Map, SOM)方法被广泛应用于地震相分析。通过识别训练好的SOM网络输出层神经元与地震相类别之间的关系，可以有效确定待测试样本与地震相类别之间的关系。SOM网络输出层神经元与地震相类别之间的关系是通过对SOM网络输出层神经元进行K均值聚类实现的。K均值聚类方法将待分析的样本集分为K个簇，使簇内样本具有较高的相似度，而簇间样本具有较低的相似度。相似度可由每个样本到簇中心的欧式距离、曼哈顿距离或马氏距离等方式度量。大量实践表明，K均值聚类易受“噪声”和孤立点的影响，无法实现对SOM输出层的神经元的最佳聚类。由此SOM输出层神经元映射到地震相类别时可能会使地震相边界模糊甚至混淆。值得注意的是，基于SOM的地震相分析方法不是利用钻井、测井、地质等先验知识，而是仅根据经验判断类别数K。张等^[4]利用粒子群优化算法对SOM输出层的神经元进行了优化聚类，但仍未利用钻井、测井等宝贵信息。

现阶段地震相分析大多基于原始地震波形。由于大地的滤波作用，地震数据呈低频、窄带特征，因而仅仅利用原始地震波形难以描述地震反射信息中的微小变化，如小断层、阻抗差异小的岩性边界。近年来，应用于油气勘探开发的各种地震属性计算方法已经取得了重大进展^[5]。众多地球物理学家们将时频分析、瞬时属性、几何结构属性(曲率、相干等)、纹理属性、弹性参数反演等结果与SOM方法相结合，实现了地震相自动划分和识别^[6]，尤其是图像纹理分析得到广泛关注。Gao^[7]将地震纹理定义为相邻地震振幅的空间排列，并指出，对于识别地震特征，与平均振幅相比，地震纹理分析有效强调并突出地震反射振幅的空间变异性。地震纹理属性已被应用于构造解释、盐丘体识别、储层描述和预测、沉积环境分析等，也被成功地应用于地震相分析^[7-8]。胡英等^[9]联合地震纹理属性和模糊C均值聚类算法实现了地震相的自动划分。但现有方法均是从叠后地震数据中提取地震纹理特征。叠前地震属性已经成为地震属性发展的一个必然方向^[10]。Song等^[11]提出基于叠前地震纹理属性地震相分析方法，并初步应用于断层和溶蚀孔洞发育区的识别。

因此，本文提出基于叠前地震纹理属性的半监督地震相分析算法，不仅能够从叠前地震数据中挖掘更多的有关微小空间、振幅随方位角/炮检距等的变异性，也能有效利用钻井、测井、地质等宝贵的先验信息，并采用全局优化算法优化SOM从输出层到地震相类别的映射关系。本文方法在理论模型数据和塔里木盆地实际资料的应用中均取得良好的效果。

1 方法原理 1.1 叠前地震纹理特征提取

相比于单一地震波形本身的特征，纹理特征不仅可以反映叠前地震反射面元本身的波形特征，还包含多个方向的空间信息，它是地震波形、空间信息、储层特征的有机整合体。在假定图像中各像素间的空间分布关系包含了图像纹理信息的前提下，Haralick等^[12]提出了分析纹理的灰度共生矩阵(GLCM)统计法，并定义了14种纹理属性。Ulaby等^[13]通过实际资料应用发现，只有能量、对比度、同质性和相关性之间是相互独立的。能量反映纹理信息的粗细程度；对比度反映纹理特征变化的程度；同质性反映局部纹理的均匀程度；相关性则反映GLCM中行元与列元的相似程度。

随着地震采集仪器的进步和观测系统的优化，共反射点叠前地震数据局部叠加的信噪比已经能够满足储层预测、流体识别的要求。根据Song等^[11]设计的叠前纹理特征提取算法，并考虑不同空间方向上的差异性，从6个方向上计算各地震反射元的GLCM。在GLCM的基础上，计算能量、对比度、相关性和同质性四种非相关的叠前地震纹理属性。本次研究不考虑时间方向，每个地震反射元的叠前地震纹理元与GLCM统计方向分别如图 1和表 1所示。

Inline和Crossline方向提取的纹理属性分别用于增强地震反射元之间在主测线和联络测线方向的关系(如断层、地质体边界、地质体横向连续性等)；炮检距方向提取的纹理属性可突出地震反射元中与岩性和流体性质相关的振幅随炮检距变化的信息；方位角方向提取的纹理属性可增强地震反射元中与各向异性(如裂缝发育程度)相关的振幅随方位角变化信息。对于不同的叠前地震数据道集和不同的应用目的，提取叠前地震纹理特征的方式存在差异，如对于描述储层各向异性特征，叠前纹理属性特征提取的重点方向是方位角方向，并且考虑到叠前纹理属性描述的是相同反射点在不同方向上地震反射振幅的变化，宽方位角叠前道集需要采用各向异性速度完成动校正处理，即保证道集被校平。

1.2 自组织映射

自组织映射是一种无监督的竞争性学习方法，由2层网络结构组成(图 2)。

图 2中输入层的维数与每个样本向量的维数一致。输入层与输出层之间采取全连接，而输出层各神经元之间为侧抑制连接。假设训练样本X_i的维度为D，首先计算X_i与竞争层各神经元之间的距离，从中找出距离X_i最近的最佳匹配神经元(Best Matching Unit，BMU)，然后根据式(1)更新神经元及其领域内包含的神经元的权值

$ \begin{array}{l} {\mathit{\boldsymbol{W}}_s}(t + 1) = {\mathit{\boldsymbol{W}}_s}(t) + \eta (t)h\left[ {{\mathop{\rm dist}\nolimits} \left( {{\rm{BMU}}, {\mathit{\boldsymbol{v}}_s}} \right), t} \right] \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {{\mathit{\boldsymbol{X}}_i} - {\mathit{\boldsymbol{W}}_s}(t)} \right] \end{array} $

(1)

式中：v_s表示领域内第s个神经单元；W_s(t)表示第s个神经元在t时刻对应的权向量；η(t)是一个随时间增加而逐渐减小的学习率；h(·)是领域函数自变量；dist(·，·)表示求取领域内神经元间的距离。依次利用训练样本更新网络权值。训练结束后，相邻的输出层节点就可以认为关联了输入数据中的相似样本。

1.3 半监督自组织映射

数据集X={x_i^u|i=1, 2, …, N_u，x_j^l|j=1, 2, …, N_l}，标签集L={l_j|j=1, 2, …, N_l, l_j∈[1, K]}，其中x^u是无标签样本，x^l是有标签样本，N_u是无标签样本数，N_l是有标签样本数，l_j为x_j^l对应的标签，K为类别数目。SOM网络权值初始化后，从X中随机抽取训练样本x输入网络中，更新X，找到x的BMU，并根据式(1)对构建的网络进行更新。重复该过程，直到X为空集。然后确定x^l在竞争层中的BMU^l，依据L对BMU^l进行标记，获得BMU^l之间的“必连” (Must-link)与“勿连”(Cannot-link)关系，最后以获得的BMU^l之间“必连”与“勿连”关系约束竞争层神经元进行优化K-means聚类，从而实现半监督自组织映射(Semi-Supervised Self-organizing Map，SS-SOM)。

1.4 算法流程

结合少量的已知标签信息，本文研究了基于SS-SOM聚类方法。为了缓解K均值聚类方法从SOM输出层到地震相映射时相边界的混淆，利用全局优化算法和监督信息进行SOM输出层关联节点的划分。采用全局优化算法对SOM输出层的神经元进行优化聚类，从而得到一个从输出层结点到地震相类别指派的最优映射，实现流程如图 3所示，具体算法步骤如下：

(1) 根据应用目标，采取合适方式从叠前地震数据中提取叠前地震纹理属性；

(2) 对训练集(包括少量有标签和大量无标签的样本)进行SOM训练；

(3) 标记标签样本对应的SOM输出层的响应元，并根据标签信息确定它们之间的“必连”或“勿连”关系；

(4) 利用全局优化算法和获取的连接约束关系对SOM输出层进行优化K均值聚类，待优化的目标函数f为K均值聚类的类内距离之和

$ f = \sum\limits_{j = 1}^{{N_{\rm{T}}}} {\sum\limits_{i = 1}^{{N_j}} {\left\| {{\mathit{\boldsymbol{x}}_{ji}} - \mathit{\boldsymbol{x}}_{\rm{c}}^j} \right\|} } $

(2)

式中：N_T为总的类别数；N_j为第j类的元素数目；x_c^j为第j类的中心；x_ji为第j类的第i个样本。注意每一次的K均值聚类必须满足步骤(3)得到的约束关系；

(5) 重复步骤(4)，直至满足全局最优化算法的停止条件，得到最终SOM输出神经元与类别的映射关系；

(6) 将待测试样本输入步骤(5)获得SOM网络，得到它在SOM输出层的响应元，根据优化的映射关系得到待测试样本对应的地震相类别。

2 模型分析

用图 4所示地质结构模型合成包含振幅随炮检距变化信息的叠前地震数据，验证叠前地震纹理属性的优越性和半监督地震相分析方法的有效性。数据生成过程中，震源子波为主频等于30Hz的雷克子波，正演合成记录采用Aki-Richards三项AVO方程，时间采样率为1ms。图 4a所示地质结构模型包含了4大套地层，不同地层岩石的物性参数(纵、横波速度、密度)如图 4b所示。图 4中存在4种横向差异岩性体。为保证合成叠前地震道集与实际采集地震数据更一致，使用带通滤波器(10~15—50~60Hz)获得信噪比为10dB的叠前地震道集数据体，并叠加生成的叠前地震道集，获得叠后地震数据体。

首先对目标数据提取叠前地震纹理属性，再利用SOM、SS-SOM方法分别生成地震相图，结果如图 5所示。SOM和SS-SOM网络结构中输出层神经元数目均为100。

图 5表明，噪声和叠加操作使叠后地震数据包含的信息模糊，致使分类结果中的类别混淆，相带完全不可识别；基于叠前地震纹理属性的地震相分析结果则相对清晰。图 5b表明，在含噪声情况下，尽管根据叠前地震纹理属性可分析出解释所需的有价值信息，但仍存在大量的错误地震相类别指派，尤其是对于第②类规模较小岩性体(黑色箭头)空间展布的预测误差更大。图 5c表明本文方法得到的地震相图误差极小。

3 实例分析

以塔里木盆地A区的巴仕基奇克组为例分析和验证本文方法的实际应用效果。巴仕基奇克组顶如图 6绿线所示。该试验区中，巴仕基奇克组地层中包含大量断层和不同厚度的砂体。根据地震资料的主频、砂岩速度以及钻井资料，将该目标地层划分为四类地震反射模式：无砂模式(如W2、W4、W5、W7井)、薄砂模式(小于调谐厚度，如W6、W9、W11井)、厚砂模式(W3、W8井)以及断层模式。

对叠前地震道集(包含振幅随炮检距变化)提取地震纹理属性(图 7)。可见，在叠前地震同质性纹理属性中，断层呈现为低同质性特征(黑色)。

借鉴de Matos等^[3]的研究成果，本文应用于实际数据处理的SOM网络结构中输出层神经元数目为100。不同方法的分析的地震相结果见图 8。

图 8a未能清晰地描述断层，并且相带的边界不清晰(如白色箭头所示)。图 8b表明，基于叠前地震纹理特征的无监督SOM聚类方法获得的地震相图对断层的空间展布描述效果有极大的提升(如红色箭头标记处)，证实了叠前地震纹理属性可有效描述实际地震数据包含的细微信息，但是地震相边界仍然模糊不清。

在进行半监督聚类的地震相分析过程中，采用交叉验证方法优选获得最终的SS-SOM网络模型，即随机地将钻遇无砂模式、薄砂模式、厚砂模式的井均留出一口井作为验证，将其余井作为标签参与训练，从而有效地避免SS-SOM模型欠学习，处理结果见图 8c。对比图 8b与8c表明，本文方法对断层和相边界的刻画都更加清晰，并且横向分辨率也得到了一定程度的提升(如图 8c中白色椭圆所示)。结合已知钻井信息对比分析，图 8a中的井W4、W6、W8和W11的信息均与测井资料不符；图 8b中井W9和W11所处位置被划分为无砂地震反射模式，井W6所处位置被划分为厚砂岩地震反射模式；而图 8c中，各井位置的相划分结果均与测井资料基本相符，证明了本文研究方法的适用性和正确性。

4 结论

本文提出了一种基于叠前地震纹理特征的半监督地震相分析方法。结合钻井已知资料和全局优化算法，采用半监督方式对传统的无监督SOM方法进行优化和改进，实现了SOM网络输出层神经元到地震相类别的更优聚类，提高了地震相带边界识别能力，能产生更准确的地震相图。数值模型和实际资料应用结果证明，基于叠前地震纹理特征的地震相分析有助于提高地震相分类的效果，对微小相带的展示更清晰和准确。