0 引言

如何高效、准确地识别断裂系统、解释断层一直是石油地球物理勘探研究的热点与难点^[1-4]，断层解释成果数据是三维地质建模的重要数据来源，可为复杂断块油藏定量描述与开发提供论证基础^[5-8]。在地震资料解释过程中逐一命名断层是一项繁琐而艰巨的工作，解释人员为了提高工作效率往往只命名一些重要断层，或者将断层分组，每组有几条甚至几百条断层。地质建模时依据不同的断层名区分不同断面，若不同断层的断层名相同，会造成虚假断面的现象。因此，在断层建模之前需要拆分断层。

手工拆分断层工作繁琐且十分耗时，而且断层拆分精度受主观因素制约，特别是在复杂断块区域，可能会造成错误的拆分结果。由于不同工区的地质问题不同，解释精度要求也不一致，致使自动拆分断层难度极大。因此，前人针对自动拆分断层技术的研究较少。田静^[9]利用深探软件，提出基于聚类分析的自动拆分断层技术，在萨尔图油田取得了较好效果。本文提出了基于张量投票的自动拆分断层方法，详细介绍了方法原理与实现过程。渤海F工区实例证实该方法的精度在90%以上，并极大地提高了工作效率。

1 基本原理

断层解释成果(断棱数据)通常包含坐标(x、y、z)、线号(Inline number)、道号(Crossline number)和断棱编号等信息，在三维空间可以视为表示位置信息的点、线数据，但拆分断层仅仅依靠位置信息是不够的。受构造运动影响，断层之间往往存在复杂的交切关系，在断层拆分时不仅需要断层的空间位置信息，还需要断层的产状(走向、倾向、倾角等)信息进行约束。因此，需要计算断层点集方向信息。计算点集数据方向信息的方法包括主成分分析法(PCA)^[10-11]、最小二乘拟合法^[12-13]、随机采样一致性(RANSAC)^[14-15]和张量投票法^[16-18]等。张量投票法具有较强的鲁棒性，因此本文选择张量投票法计算断层点集的方向属性。

三维张量投票中的张量为一个非负定三维对称方阵，其特征值表示一个点位于某种几何结构的可能性(即显著性)，对应的特征向量则代表该几何结构的方向。三维张量T在几何上描述为一个椭球体(图 1a)，对其特征值分解可了解几何意义

$ \begin{gathered} \boldsymbol{T}= {\lambda _1}{{\boldsymbol{e}}_1}{\boldsymbol{e}}_1^{\text{T}} + {\lambda _2}{{\boldsymbol{e}}_2}{\boldsymbol{e}}_2^{\text{T}} + {\lambda _3}{{\boldsymbol{e}}_3}{\boldsymbol{e}}_3^{\text{T}} = \hfill \\ \;\;\;\left( {{\lambda _1} - {\lambda _2}} \right){{\boldsymbol{e}}_1}{\boldsymbol{e}}_1^{\text{T}} + \left( {{\lambda _2} - {\lambda _3}} \right)\left( {{{\boldsymbol{e}}_1}{\boldsymbol{e}}_1^{\text{T}} + {{\boldsymbol{e}}_2}{\boldsymbol{e}}_2^{\text{T}}} \right) + \hfill \\ \;\;\;{\lambda _3}\left( {{{\boldsymbol{e}}_1}{\boldsymbol{e}}_1^{\text{T}} + {{\boldsymbol{e}}_2}{\boldsymbol{e}}_2^{\text{T}} + {{\boldsymbol{e}}_3}{\boldsymbol{e}}_3^{\text{T}}} \right) \hfill \\ \end{gathered} $

(1)

式中λ₁＞λ₂＞λ₃为特征向量e₁、e₂、e₃对应的特征值。式中第一项对应三维棒张量，代表一个法向量为e₁的平面。第二项对应一个圆盘状椭圆体，被命名为板张量。第三项对应一个三维球张量，代表一个没有主方向的结构(图 1b)。

在构造点的初始张量表示后，每一个点将在邻域范围内产生张量场。从实现层面上看，“投票”过程就是空间邻近点产生的张量场互相叠加的过程。所谓张量场叠加，表示投票点在邻域产生的张量乘以投票点的显著性并与接收点张量矩阵直接相加。通过张量投票，接收点的张量经叠加，使特征值和特征向量方向发生变化，意味着接收点的方向和显著性不同于初始值。这一过程实现了邻域内结构信息的传递。

由于存在不同形式的张量，计算张量投票的张量场也分为不同形式。Wu等^[19]将棒张量场、板张量场和球张量场统一，推导了闭合形式的张量解。在空间数据中，在点x_j的张量K_j投票到点x_i的张量为

$ {{\boldsymbol{S}}_{ij}} = {c_{ij}}{{\rm{\boldsymbol{R}}}_{ij}}{{\boldsymbol{K}}_j}{\rm{\boldsymbol{R}}}{{\rm{'}}_{ij}} $

(2)

式中：K_j为二阶对称张量，由张量类型决定；R_ij=I-2r_ijr_ij^T，I为单位矩阵，r_ij为从点x_j指向点x_i的单位向量；$ {\rm{\boldsymbol{R}}}{{\rm{'}}_{ij}} = \left( {{\boldsymbol{I}} - \frac{1}{2}{{\boldsymbol{r}}_{ij}}{\boldsymbol{r}}_{ij}^{\text{T}}} \right){{\rm{\boldsymbol{R}}}_{ij}} $；而

$ {c_{ij}} = \exp \left( { - \frac{{{{\left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_j}} \right\|}_2}}}{{{\sigma _{\rm{d}}}}}} \right) $

(3)

式中σ_d为尺度参数。经过张量投票后，在点x_i的张量通过叠加所有投票到该点的张量场得到

$ {{\boldsymbol{K}}_i} = \sum\limits_{{{\boldsymbol{x}}_j}} {{{\boldsymbol{S}}_{ij}}} $

(4)

式中x_j为在点x_i邻域范围内的点。

由于原始断层点集没有方向信息，在张量表示时所有点的初始张量都被定义为单位矩阵，即K_j=I。然后利用式(2)、式(4)计算、更新张量场。当所有点完成张量计算、更新之后，对K_i使用特征值分解，得到按降序排列的特征值λ_i1、λ_i2和λ_i3以及对应的特征向量e_i1、e_i2和e_i3。其中e_i1=(e_1x, e_1y, e_1z)为x_i点所在断层面的法矢量。根据e_i1可以计算断层点集中每个点对应的瞬时倾角θ和瞬时倾向Ψ

$ \left\{ {\begin{array}{*{20}{l}} {\theta = \arccos {e_{1z}}}\\ {{\mathit{\Psi}} = {\rm{arctan}}\frac{{{e_{1x}}}}{{{e_{1y}}}}} \end{array}} \right. $

(5)

利用张量投票法计算断层点集方向的流程(图 2)为：①将所有散点编码，并以球张量的形式表现；②根据σ_d计算所有散点的球张量场；③球张量投票；④对投票后的张量场进行特征值分解，得到所有散点的法向量；⑤根据法向量计算所有散点的θ和Ψ。在利用张量投票法计算断层点集方向属性时，唯一需要控制的参数为σ_d。

2 断层产状参数分析

为了验证方法的效果，并分析σ_d对方向属性计算的影响，运用模拟交叉断层数据和实际工区数据分别测试。

首先，根据存在交切关系的断层，模拟生成了相交断层点集模型(图 3a)，分别设定不同的σ_d得到张量投票法计算结果(图 3b~图 3f)。分析计算结果可知：当σ_d取值过小时(1、5m)，离散点之间无法建立投票通道，因此无法计算方向信息(图 3b、图 3c)；当σ_d取值合适时(10、20m)，主要由同一个平面的离散点相互建立投票和通信通道，因此较准确地计算了方向信息(图 3d、图 3e)；当σ_d取值过大时(100m)，同一个平面和不同平面之间的离散点均能相互建立投票和通信通道，导致计算的方向信息偏差较大(图 3f)。此外，实际断层点集的三维空间分布与模拟数据不同，纵向分布往往比横向稀疏很多，同一断棱上相邻点的纵向距离可达600m以上，而同一断层的相邻断棱的水平距离与解释的线间隔有关，一般为125~250m。测试工区相邻断棱的平均水平距离为125m，为了保证断层点集的纵向距离、水平距离的均匀性，首先对断层点集纵向插值，之后抽取部分断层点集按照不同的σ_d计算方向信息(图 4)。可见：当σ_d为50m(图 4a)时，利用同一断层相邻线的断棱无法建立投票和通信通道，此时计算误差较大；当σ_d为300m和500m时计算结果较准确(图 4b、图 4c)；当σ_d为2000m(图 4d)时，部分断层的边界点受其他断层的影响，导致较大角度误差。理论数据和实际数据测试结果表明，当σ_d设定为横向断层点平均距离的2~4倍时，能较好地计算断点产状，为后续断层拆分提供依据。

3 断层拆分 3.1 初步拆分

通过计算得到断点的方向信息后，可以利用的信息包括空间坐标(x，y，z)、断棱编号(Segment ID，可判断断点是否属于同一条断棱)、Inline号和Crossline号、断层方向信息(瞬时倾角、瞬时倾向和走向)。根据这些信息，首先初步拆分断层，具体拆分流程为：①从断层点集中随机选取同一条断棱的点集作为初始断面模型；②从余集选取小于预设的距离阈值(同一断层相邻断棱距离的2~4倍)、倾角阈值和倾向阈值的点，构成新的断面模型；③修正新的断面模型的平均法矢量，重复步骤②，直至找不到新的样本点为止；④对找到的断层模型标号，并从断层点集删除该断棱数据，然后重复步骤①；⑤算法结束。

3.2 精细拆分

在初步拆分过程中将具有相近空间位置和相似断层产状的数据分为一类，形成一条断层。这种方法在构造相对简单的工区可以有效拆分断层。但是遇到空间上非常接近且具有相近产状的多条断层时，使用正常的距离和方向阈值不能初步拆分断层(图 5a)；如果调整距离和方向阈值，则可能造成过度拆分。因此需要结合地质认识对初步拆分结果质控，若不满足需求，则对所有初步拆分结果精细拆分。

在初步拆分的基础上，融合地震、地质解释的先验知识：①同一断面在每条Inline测线通常只有一条断棱；②同一断层的局部产状变化不大，即相邻线的断棱应具有相近的倾角，并且局部走向与断层总体走向相差不大；③同一断层在Inline和Crossline方向的断棱闭合差较小。经过精细拆分(图 5b)，将相邻且产状相近的断层(图 5a)成功拆分成为两条断层。

4 实际数据测试

在实际断层拆分时，相邻断层距离越小、产状越相近，拆分难度越大。因此，对于同一工区，先对一部分断层分组命名可大大减小拆分难度，后期可以依次人工拆分(只是耗时较长)。否则，当构造复杂时，手工拆分断层的难度非常大。本文选取构造十分复杂的F区，测试在近似极限情况下的断层自动拆分效果。F区位于渤海湾盆地西部，面积约为800km²。受多期次构造运动影响^[20]，发育三组左阶排列走滑断层(NE向、NWW向和近EW)，断层之间距离很近且产状相近，存在“Y”字型相交断层，剖面上呈花状构造^[21-25]。为了测试效果，将所有原始断层数据命名为一个名字(图 6a中所有断棱颜色相同)，在此基础上拆分断层。可见：①与原始数据(图 6a)相比，初步拆分结果(图 6b)删除了部分小规模孤立断棱，能准确拆分大部分断层，但部分断层拆分不彻底(图 7上左、图 7下左)。主要原因是这些断层空间接近、产状相近，甚至出现部分相连的情况，因此在相同参数情况下，仅利用初步拆分算法很难一次性完成准确拆分。②精细拆分结果基于初步拆分结果，将初步拆分不彻底的断层进一步拆分(图 6c，图 7上右、图 7下右)，断层位置更准确、合理，更有利于建模，精细拆分结果的准确率大于90%。

5 结论与建议

本文提出了一种自动拆分断层方法，该方法的核心是利用张量投票算法计算原始断层点集的方向信息，然后结合断层文件中的空间位置和断棱编号信息，融合现有地质认识，实现断层的自动拆分。精细拆分融合了地震、地质解释的先验知识，通过走向差阈值和闭合差阈值的约束，进一步拆分初步拆分结果。对于走向差阈值和闭合差阈值的选取，则需要建模人员以工区地质背景为参考。断层自动拆分可以极大地提高建模效率，以F区为例，传统的手动拆分耗时在1周以上，而本文提出的自动拆分算法在准确率大于90%的前提下，可在30分钟内完成拆分，因此可以作为一种有效的建模辅助手段，使建模人员摆脱繁琐的手动拆分工作，同时为地震资料解释人员提供一种高效、快速拆分、命名断层的方法。下一步的研究方向是减少阈值的硬约束，进一步提高自动拆分算法的智能性。