基于波形特征向量的凝聚层次聚类地震相分析

刘仕友，宋炜，应明雄，孙万元，汪锐

(1.中海石油(中国)有限公司 湛江分公司，广东 湛江 524057； 2.中国石油大学(北京) 地球物理学院, 北京 102249)

0 引言

Sloss等[1]认为地震沉积相表达的是沉积物在地震记录上的响应特征；Vail等[2]指出地震相分析是利用地震数据来解释地层的沉积环境和岩相过程的一门技术。曾洪流等[3-8]提出了有关地层切片的方法、理论模型研究和实际地震资料解释与应用效果，并在“The Leading Edge”上首次使用了“地震沉积学”这个名称[9]。自20世纪90年代起，大量研究证实地震地貌学是沉积成像研究的有力工具，进而衍生出地震沉积学。目前，较为先进的地震相划分方法是基于机器学习和深度学习的地震相定量分析方法，主要包括两种：有监督聚类和无监督聚类分析[10-11]。有监督聚类方法需要根据有限的钻井资料和地震数据建立训练分类器，然后使用该训练器对无井区域的地震数据进行预测。但是由于训练样本较少，地震波形或者地震属性与已知含气性很难建立较好的关系，这样预测的精度会降低。刘爱群等在分频属性分析确定砂体边界的基础上，结合已知油气田的波形特征，开展聚类分析进一步确定有利储层[12]。白博等则是在计算能够突出贝壳灰岩内部反射特征的伪阻抗体的基础上，开展伪阻抗波形聚类分析预测贝克灰岩展布[13]。李辉等也基于自组织神网络波形聚类的方法来识别出主河道、河道侧翼等不同叠合模式砂体边界[14]。徐海等提出了高精度层序地层格架约束下的基于波形相对变化的波形-微相定量表征综合解释技术，首先是利用主成分分析法PCA(principal component analysis)降维，然后采用有监督和无监督自组织神经网络方法SOM(self-organizing map)实现波形聚类分析[15]。张艳等提出了基于SOM和HSV染色技术联合开展致密砂岩储层地震相分析方法[16]。学者们主要采用较为传统的聚类分析方法来进行研究工作，取得了一定的应用效果，但并没有将地震沉积学原理的地震相分析思想有机地结合到实际应用中。在汲取前人研究基础上，我们引入了一种新的基于波形聚类的划分地震相的分析方法。首先设计了一个三维薄层的层间地质模型，通过物理模拟得到三维地震记录；然后，利用常规的基于地震沉积学原理的地层切片技术开展目标层地震相分析；最后，采用基于无监督机器学习的波形特征向量凝聚层次聚类开展地震相分析工作。该方法与传统的地层切片地震相分析方法相比较，不仅考虑了地震信号的振幅属性特征，而且通过波形特征的变化，综合考虑了地震信号的振幅、相位和频率属性特征。物理模型试验和实际地震资料应用表明，该方法能够更准确、更合理地描述地震相的空间分布。

1 基本原理

1.1 凝聚层次聚类分析

原型聚类假设聚类结构能够通过一组原型刻画，先对原型进行初始化，然后对原型进行迭代更新求解。层次聚类(hierarchical clustering)是一种基于原型的聚类算法，试图在不同层次对数据集进行划分，从而形成树形的聚类结构。数据集的划分可采用“自底向上”的聚合策略，也可以采用“自顶向下”的分拆策略。层次聚类算法的优势是以通过绘制树状图(dendrogram)使用可视化的方式来解释聚类结果。层次聚类可以分为凝聚(agglomerative)层次聚类和分裂(divisive)层次聚类。分裂层次聚类采用的就是“自顶而下”的思想，先将所有的样本都看作是同一个类簇，然后通过迭代将簇划分为更小的簇。凝聚层次聚类采用的是“自底向上”的思想，先将每一个样本都看成是一个不同的类簇，通过重复将最近的一对类簇进行合并。在凝聚层次聚类中，判定类簇间距离的两个标准方法就是单连接(single linkage)和全连接(complete linkage)。单连接是计算每一对簇中最相似两个样本的距离，并合并距离最近的两个样本所属类簇。全连接是通过比较找到分布于两个类簇中最不相似的样本(距离最远)，从而来完成类簇的合并。凝聚层次聚类除了通过单连接和全连接来判断两个簇之间的距离外，还可以通过平均连接(average linkage)，平均连接定义簇的邻近度为取自两个不同簇的所有点对邻近度的平均值。图1是不同连接条件下随机生成的样本点分成9类的效果的差异，图中不同颜色代表不同类别，其中数据连接类型分为：图1a是单连接，图1b是平均连接，图1c为全连接，可见连接类型不同，聚类结果也不同。在实际计算中，一般都是将距离最近的两个簇进行合并，对于两个点，距离比较好计算，欧氏距离或者闵可夫斯基距离都可以。但是两个簇(簇中不仅包含一个点)的距离算法有多种，其步骤基本相同，差别在于聚类间距的定义不同。计算类簇距离有如下 4个广泛采用的度量方法，其中|P-P′|是两个对象P和P′之间的距离，mi是簇Ci的均值，mj是簇Cj的均值，而ni是簇Ci中对象的数目，nj是簇Cj中对象的数目。这些度量有称连接度量(linkage measure):

最小距离：

(1)

最大距离：

(2)

均值距离：

distmean(Ci,Cj)=|mi-mj| ,

(3)

平均距离：

(4)

a—单连接;b—平均连接;c—全连接a—single connection;b—average connection;c—full connection图1 相同数据集不同连接类型的分类结果Fig.1 Classification results of different connection types in the same data set

一些经典的聚类算法如k-means、密度聚类等算法都要先随机指定k个类别中心点，类簇用这些中心点逐步迭代找到最终的中心点，如果随机选取的点不具有代表性，会导致聚类效果比较差。凝聚层次聚类算法不需要提前指定随机的中心点，就能够有效地对初始样本进行聚类，但也需要指定聚类类别数k或者终止条件。常规的聚类算法如k-means等在样本维度较高，样本数较大时，会导致聚类分析时间过长等问题，而凝聚层次聚类则适用于较大的数据集的聚类分析，尤其是在有连接约束(connectivity constraint)的时候。图2是随机生成的空间分布样本点，利用本算法聚成6类，图中不同颜色代表不同类别，左边是无连接约束(without connectivity constraint)，可以看到有些青色的簇横跨了两片，如红色箭头指；黑色的簇横跨了两片，如黑色箭头指。右边是有连接约束(with connectivity constraint)，类簇是沿着弯曲的平面分布的，这种结果更合理。对于有连接约束，在聚类样本点很大的情况下，对于每个点只需要考虑和它相邻的点，而不是考虑所有的点，因此可以加快计算速度。在地震波形聚类划相时，考虑相带的连续性和变化，选择有连接约束会更好。

a—无连接约束;b—有连接约束a—connectionless constraint;b—connectivity constraints图2 相同数据集不同连接约束Fig.2 Different connection constraints for the same data set

1.2 凝聚层次聚类算法过程

输入：沿地层切片提取地震波形特征向量构成样本集{x1,x2,…,xn}，选择类簇距离度量函数dist，给定类簇个数k。

算法过程：

①每个样本属于单独的一类，即Ci={xi};

②计算样本与样本之间的距离，用矩阵M表示，其中Mij=dist(Ci,Cj);

（1）经检测可知，120例患者的冠状动脉支共390支，其中324支的图像质量为1级，56支的图像质量为2级，15支的图像质量为3级。

③当前类簇数为q=n，开始簇合并聚类;

④循环while(q>k);

a.找出距离最近的两个簇Ci和Cj，令Ci=Ci∪Cj

b.将Cj后面的簇的编号向前移动一位

c.删除矩阵Mij=dist(Ci,Cj)的第i行和j列

d.更新Mij=dist(Ci,Cj)矩阵(主要更新i簇与其他簇的距离),j=1:q-1

e.更新当前类簇个数q=q-1

⑤输出生成的簇{C1，C2，…，Ck} 。

2 物理模型设计及凝聚层次聚类方法应用效果

图3 物理模型采集尺寸及参数示意图Fig.3 Acquisition size and parameter sketch of physical model

图4 物理模型6层砂体空间分布示意图Fig.4 Spatial distribution diagram of six Layers of sand body in physical model

图5 砂体叠置模型垂直切面示意Fig.5 Diagram of vertical section of sand overlay model

地震沉积学义时间平行的参考同相轴代表了地质层面或沉积单元，控制着地震反射同相轴的反射能量。图6是物理模型正演模拟三维数据体中一条地震剖面，图中的彩色层位线是根据4个控制层划分的地层切片。常规的基于地震沉积学原理的地层切片均方根振幅属性在沉积相划分时可以有效提高地震数据的横向分辨率，但是由于均方根振幅属性没有综合考虑地震波形中频率和相位信息的影响，因此，其划分地震相的精度是有限的，且容易受噪声的影响。相比而言，采用无监督机器学习地震波形凝聚层次聚类方法，由于利用波形信息进行聚类分析，综合考虑了振幅、频率、相位信息，能有效克服单纯基于均方根振幅属性划分地震相的缺点，提高划相的精度。

机器学习算法一般需要假设输入数据是正态分布的(即满足具有零均值和单位方差的高斯分布)，因此需要对输入数据进行标准化，使其满足正态分布。在实际资料处理中，首先基于地震沉积学原理制作地层切片，沿地层切片在一定时间窗范围内提取波形；然后，采用机器学习算法包Scikit-learn中的StandardScalar函数对输入的地震波形数据进行标准化；最后，通过波形特征向量凝聚层次聚类将不同的波形分为不同的类簇。图7是沿图6中黄色箭头所指的地层切片提取的200个波形曲线，其波形维度为31。

图6 物理模型偏移剖面及基于地震沉积学原理解释的层位Fig.6 Physical model migration profile and horizon interpretation based on seismic sedimentology principle

图7 凝聚层次聚类输入的200个波形曲线Fig.7 200 waveform curves inputted by AHC

图8a是第一层“蛇形”砂体平面形态，并将砂岩顶底叠置显示，其中黄色部分是砂岩底的范围，白色部分是砂岩的顶。在实际制作中，由于模型误差和环氧树脂中添加滑石粉的过程中无法保证能够混合均匀，因此会产生速度的变化。图8b是沿地层切片提取的均方根振幅属性，可以较好地描述砂体的展布和形态，但无法分辨出砂体顶底。图8c是波形凝聚层次聚类结果，聚类类别参数K为7, 波形向量31维。图8d是波形凝聚层次聚类结果，聚类类别参数K为7, 波形向量37维。波形凝聚层次聚类分析结果表明，本方法能够很好地描述砂体的空间展布及顶底边界的范围，波形的维度会明显影响分类的结果，因此，确定合理的聚类分析波形维度关键参数，可获得较好的聚类结果。

物理模型中其他四层砂泥岩的空间展布及形态如图9所示，其中图9a是第二层“肠状及指状”砂体平面形态；图9b是第三层“肠状”砂体和“椭圆形”泥岩平面形态展布；图9c是第四层“哑铃状”砂体和“指状”泥岩平面形态展布；图9d是第五层“菱形”砂体和“点状”泥岩平面形态展布。图10是沿地层提取的均方根振幅属性切片：图10a“肠状”砂体和“指状”砂体的平面展布清晰可见，但指状砂体正下方的大片薄砂体无法准确描述；图10b“肠状”砂体和“椭圆形”泥岩的空间展布形态可以清楚地划分；图10c“哑铃状”砂体和“指状”泥岩刻划的比较清楚，只是在“哑铃状”砂体的正上方，有一部分蓝色的振幅属性，模糊了“哑铃状”砂体的边界；图10d大一点的点状泥岩可以清楚地刻划出来，三个菱形砂体比较清晰地显现，但下方的边界不清晰，正下方的“点状”泥岩未显现，菱形砂岩上方偏南的“点状”泥岩也未显现。

a—第一层“蛇形”砂体平面形态;b—沿地层切片提取的均方根振幅属性;c—波形凝聚层次聚类结果(K=7，波形向量31维);d—波形凝聚层次聚类结果(K=7，波形向量37维)a—the plane shape of the first layer of "serpentine" sand body;b—RMS amplitude attributes extracted along stratum slices;c—waveform AHC(K=7,waveform vector 31 dimensions);d-waveform AHC(K=7,waveform vector 37 dimensions)图8 物理模型第一层砂体平面形态、均方根振幅属性及波形凝聚层次聚类结果Fig.8 Sand body planar morphology,RMS amplitude attributes and waveform AHC in the first layer of physical model

a—“肠状及指状”砂体平面形态;b—“肠状”砂体和“椭圆形”泥岩平面形态;c—“哑铃状”砂体和“指状”泥岩平面形态展布;d—“菱形”砂体和“点状”泥岩平面形态展布a—planar morphology of "intestinal and finger" sand bodies;b—planar distribution of "intestinal" sand body and "elliptical" mudstone;c—horizontal morphological distribution of dumbbell sand body and finger mudstone;d—plane morphological distribution of "rhombic" sand body and "point" mudstone图9 物理模型第二至第五层砂、泥岩空间展布及形态Fig.9 Spatial distribution and shape of sand and mudstone in the second to fifth layers of physical model

图11a是第二层基于波形凝聚层次聚类结果(K=7，波形向量13维)的“肠状及指状”砂体,图中“肠状”和“指状”砂体的平面展布也清晰可见；指状砂体正下方大片薄砂体也有所反映；图11b是第三层基于波形凝聚层次聚类分析结果(K=7，波形向量13维)的“肠状”砂体和“椭圆形”泥岩。从聚类分析图上看，“肠状”砂体和“椭圆形”泥岩平面形态展布特征比均方根振幅属性更加精细。遗憾的是,无论均方根振幅属性还是聚类分析属性，都无法区分“肠状”砂岩和“椭圆形”泥岩的岩性。图11c是第四层基于波形凝聚层次聚类分析结果(K=7，波形向量21维)的“哑铃状”砂体和“指状”泥岩，“哑铃状”砂体和“指状”泥岩边界都能清晰刻划，同样的问题是这两种方法都不能有效区分砂岩和泥岩。图11d是第五层基于波形凝聚层次聚类分析结果(K=7，波形向量21维)的“菱形”砂体和“点状”泥岩,从图中可以看出，三个菱形砂岩的边界可以清晰地刻划，右边的边界也能刻划的比较清楚，正下方 “点状”泥岩也清晰可见，但是菱形砂岩上方偏右的“点状”泥岩同样也没有显示。从上述分析可以看到尽管模型的层非常薄，但通过地层切片得到的地震相具有非常高的纵横向分辨率。但由于模型在浇铸过程中并不能保证砂体均匀成型，所以，砂体的背景介质并不是均匀的，因此，除了砂体，在其他部位也存在相带变化。基于均方根振幅属性地震相划分中，相带变化的边界并不是特别清晰，没能显示出某些局部相带变化。与基于均方根振幅属性地层切片得到的地震相相比，波形凝聚层次聚类的方法在地震相描述方面更为详细，地震相带的变化划分更精确，因此，基于地震沉积学原理再结合本文提出的无监督机器学习的波形凝聚层次聚类方法，可进一步提高地震相划分的精度。

a—均方根振幅属性“肠状”砂体和“指状”砂体的平面展布;b—均方根振幅属性“肠状”砂体和“椭圆形”泥岩的平面展布;c—均方根振幅属性“哑铃状”砂体和“指状”泥岩的平面展布;d—均方根振幅属性“菱形”砂体和“点状”泥岩的平面展布a—RMS amplitude attribute plane distribution of "intestinal" sand body and "finger" sand body;b—RMS amplitude attribute plane distribution of "intestinal" sand body and "elliptical" mudstone;c—RMS amplitude attribute "dumbbell" sandbody and "finger" mudstone plane distribution;d—RMS amplitude attribute "rhombic" sand body and "point" mudstone plane distribution图10 各层地层切片均方根振幅属性地震相划分Fig.10 Seismic facies division of RMS amplitude attribute of stratum slices

a—AHC属性“肠状”砂体和“指状”砂体的平面展布;b—AHC属性“肠状”砂体和“椭圆形”泥岩的平面展布;c—AHC“哑铃状”砂体和“指状”泥岩的平面展布;d-AHC属性“菱形”砂体和“点状”泥岩的平面展布a—AHC attribute plane distribution of "intestinal" sand body and "finger" sand body;b—AHC attribute plane distribution of "intestinal" sand body and "elliptical" mudstone;c—AHC attribute "dumbbell" sandbody and "finger" mudstone plane distribution;d—AHC attribute "rhombic" sand body and "point" mudstone plane distribution图11 各层凝聚层次聚类地震相分布Fig.11 Seismic facies division of AHC attribute of stratum slices

3 基于AHC的储层类别划分海洋地震数据应用实例分析

实际数据来自中国海洋石油总公司(CNOOC)位于南海西部的湛江分公司琼东南盆地深水区。目标储层是具有强反射特征的砂岩储层，在地震剖面上显示为亮点反射。研究过程中，首先利用解释的特征反射格架层，T50、T60、T70、T100为基础，利用地震沉积学原理，细分成相应的地层切片。图12是目标区基于地震沉积学原理解释的层位，图13是目标层T70等T0图，从图中可以看出构造高、低部位的空间展布。图14是叠前时间偏移地震数据沿目标层T70的均方根(RMS)振幅属性地层切片。从图中可见在构造高部位和斜坡带都有强的振幅分布，其中A部位是已经钻井证实的具有高含气饱和度的砂岩储层。在实际生产中基于均方根振幅属性，后续又确定了B、C这两个和A部位具有相同均方根属性特征的目标，都表现为较强的均方根振幅属性异常。由于强反射能量的特征差异影响因素较多，因此单凭均方根属性特征划分储层类型还是有一定风险的。在图14所示的均方根振幅属性图中，除了A、B、C三个区域表现为强振幅属性外，在相邻区域还有D、E、F等多个强振幅区域，这就为甄别有效储层造成了困难。为了更直接地划分不同反射特征的储层类别差异，采用本文提出的波形凝聚层次聚类分析方法，对沿T70地层切片提取的波形进行聚类分析(聚类类别参数K=7，波形向量21维)。

图12 目标区基于地震沉积学原理层位解释Fig.12 Horizon interpretation based on seismic sedimentology principle in target area

图13 T70层等T0图Fig.13 T70 layer isotropic T0 graph

图14 沿T70提取的叠前时间偏移数据体地层切片均方根振幅属性Fig.14 RMS amplitude attribute of stratigraphic slices extracted along T70 from prestack time migration data

图15是叠前时间偏移体目标层地震数据波形凝聚层次聚类属性，该聚类分析的结果将A、B、C三个区块分为三种不同的类型，A、B两区域的类别更相近，因为已经通过实钻确定了A区域是含气砂岩油气藏，因此，可以判定B区域将是下一步开展钻探的有利目标区。而C区域类别则是完全不同于A区域的类别，因此不能确定是下一步钻井的有利目标。此外，其他的强振幅区域，如图15所示的D、E、F区，通过聚类分析划分为不同的类型，而这些区域所划分的类型和A区域也不同， 且处于构造的低部位，因此，可以初步判定不是下一步钻探的有利目标。本方法有效地解决了利用均方根属性无法解决的储层类别差异问题，但要想进一步确定储层特征，还需要其他的技术方法辅助来开展分析研究。

图15 沿T70地层切片从叠前时间偏移数据体提取地震波形获得的凝聚层次聚类属性Fig.15 AHC attributes obtained by extracting seismic waveforms from prestack time migration data along T70 stratigraphic slices

4 结论

基于凝聚层次聚类的波形聚类地震相划分的优点在于：

1)具有更好的抗噪性；

2)物理模型试验表明，在合理选择参考等时面的情况下，基于地震沉积学原理结合本文引入的波形特征向量凝聚层次聚类分析技术可显著提高地层切片的横向分辨率；

3)与均方根振幅属性比较，凝聚层次聚类属性可提高地震相划分的可靠性，增强地震相类别划分的精度；

4)具有较高的横向分辨率，是储层相变分析和沉积相研究的有效辅助工具。

实际应用中，需要注意以下几点：

1)参考等时面的解释和地层切片的划分至关重要，将直接影响聚类分析结果的精度和可靠性；

2)用来聚类分析的波形向量维度是和储层的厚薄有关，需要针对储层的特征，沿地层切片选择时窗，合理选择波形样点数；

3)聚类类别数的设定也对最终聚类结果有一定的影响，需要结合实际情况进行测试分析。