一种面向三维层面网格生成的层序分析方法

摘要

本发明公开了一种面向三维层面网格生成的层序分析方法，用于地质层面建模和地下结构分析。其特征在于，包括：对离散点构成的面片进行扫描；根据扫描的结果建立目标函数；求取使目标函数取最小值的解，并用该解作为面片的序值；利用该序值为面片分类，并将分好类的面片作为层面网格生成的数据。本发明可实现同相轴自动追踪结果的自动分类，从而完成自动化的解释，摆脱繁重的人工分类拾取，有效提高从地震数据体到层面网格建模的解释过程的效率。

说明书

技术领域

本发明涉及利用同相轴追踪的结果进行地震数据的自动解释，并将解释的结果用于三维层面网格生成。该解释方法基于层序分析。

背景技术

石油资源是当今世界经济发展依赖的主要能源形式，同时也为许多化工行业提供原材料。可以说石油是一个国家经济发展的命脉。做为从地下获取石油资源的唯一手段，油气勘探在整个石油工业中处于基础地位。由于我国地质条件的复杂和石油资源的匮乏，我国的油气勘探面临着更加困难和艰巨的局面。

现代油气勘探一般分为如下几个环节：地质调查，野外地球物理勘探作业，地震成像，地震数据解释，基于油藏数值模拟的油气开采规划和管理。

成功的油气开采的第一步是找到油气的储藏位置，这需要我们对地下的地质构造有着精确的认识，而这一点目前是通过高精度的地震成像来完成。当前复杂地质构造成像效果最好的方法是叠前深度偏移成像。叠前深度偏移成像质量的优劣取决于我们建立的地下速度模型的好坏。由于速度与深度成像的耦合关系，得到一个精确的速度模型通常需要多轮的迭代分析，这个过程称为速度分析。在速度分析中，需要有一个地下结构的地质模型。例如在层析速度分析中，地层模型给射线追踪提供反射点和地层倾角等信息，在一些复杂的速度分析方法中，地层模型还提供约束条件。

找到油藏后，还需要根据油气运移等理论，在地下构造模型的基础上进行油藏数值模拟，根据模拟的结果设计油井及开采方式；已经投入生产的油田则利用油藏数值模拟安排生产计划。

在以上两个环节中，地质模型起到了关键的作用，反映了我们对地下结构的认知，是地球物理勘探结果的最直接体现。建立地质模型所需要的输入数据是地层离散点，其来源是地震成像数据体的解释结果。地震数据解释是指根据成像结果，分析地下结构，将地下结构信息提取出来。例如将每一个成像剖面上的代表地层分界面的强信号拾取出来，最后把所有剖面的结果综合在一起，将同属一个反射面的合并，就构成了层面离散点集。当前这一过程完全依赖于人工操作。由于数据体非常庞大，所以非常耗费人力，解释一次的周期非常长，且无法做到每一个剖面都解释，只能抽取其中的部分进行作业。所以在各个环节中，地震数据解释的次数不多，模型更新的速度小于速度更新的次数。由于速度的变化会导致成像结果的变化，理想情况下，每一次生成新的成像，都应该根据新的成像进行新的解释，生成新的地质模型，然后再进行下一轮的迭代分析。由于繁琐的人工解释，不可能达到这一要求。

为了能够将解释过程自动化，加快解释的速度，避免繁琐枯燥的作业，人们开发出了许多同相轴自动追踪算法，试图代替人工解释。同相轴的追踪可以看作是在数据场中，以某一出发点作为基准(该点称为种子点)，不断在四周寻找和种子点相似的点，扩展该集合，直到不满足连续性条件(通常为在一定邻域内找不到满足相似度的点)，然后再从新的种子点开始创建新的集合。然而，由于地下结构的复杂性和算法自身的局限性，目前同相轴自动追踪尚无法代替人工拾取，其主要原因在于自动算法生成的结果不连续，中断很常见，形成了大量零碎的面片，失去了彼此间的联系，依然需要人工分析然后将其中同属一层的面片合并，才能用于建模，所以只能作为辅助工具。

发明内容

本发明要解决的问题是：将解释的过程自动化，减少人工作业时间。由于地质建模的自动化程度已经很高，所以解释的自动化会缩短从地震数据体到地质模型的工作时间，从而提高地震数据处理的效率。

本发明从全局出发，以同相轴自动追踪的结果为基础，通过对面片集合进行层序分析，恢复追踪结果的空间拓扑信息，并利用该结果进行层面网格生成，从而实现地震数据体到层面网格的自动化。

本发明解决其技术问题所采用的技术方案是：提供了一个系统，利用同相轴自动追踪的现成结果，将地震数据的解释从局部问题变为全局问题，为同相轴自动追踪出的面片或线段赋予一个数量特征，该数量特征刻画了面片或线段在空间分布的竖直位置特征，使得同属一层的面片或线段在该特征上数值上比较接近，从而可以将面片或线段进行分类合并，直接用于层面网格建模。其特征在于，包括：

步骤A，对输入的同相轴追踪结果进行竖直方向的扫描，生成序关系计数矩阵和累积间距矩阵；

步骤B，在序关系计数矩阵和累积间距矩阵的基础上建立目标函数；

步骤C，求解使目标函数达到最小值的解，该解向量的分量为对应面片的序值；

步骤D，根据序值对面片进行分类；

上述所述面向三维层面网格生成的层序分析方法，其特征在于，所述步骤A中输入的同相轴追踪结果为对地震数据体进行同相轴追踪后得到的面片集合或线集合(三维数据体对应的是分布在三维空间的面片集，二维数据体对应的是面片集与竖直剖面的交线集)，其中面片或线是由追踪得到离散点组成，其表示的是较为连续的地层面。同相轴追踪是指人工或程序自动寻找地震数据体中在一定空间范围内连续的表示同一地层的信号，权利要求1所述的层序分析方法对同相轴追踪方法不做限定，需要的是其输出结果。

上述所述面向三维层面网格生成的层序分析方法，其特征在于，所述步骤A进一步包括：步骤A1，用竖直方向的扫描线扫描面片集合或线集，即用垂线穿过面片或线。步骤A2，每一条扫描线穿过后记录如下信息：不同的面片或线在竖直方向上被穿过的先后次序，并在扫描的过程中记录这些次序出现的次数，这些次数构成序关系计数矩阵；计算被穿过的面片或线在竖直方向上的间距，并在扫描的过程中记录间距的累积和，该累加和构成累积间距矩阵，并利用该累积和和次序出现的次数得到平均化后的间距值。

上述所述面向三维层面网格生成的层序分析方法，其特征在于，所述步骤A中，可以根据需要添加预解释的完整的控制层，扫描的时候也包括控制层在内，从而达到控制排序的目的，提高排序的精度。

上述所述面向三维层面网格生成的层序分析方法，其特征在于，所述步骤A中，除使用两面片的竖直距离外，还可以使用该面片处于上下两个面片之间的相对距离，或定义其它的距离。

上述所述面向三维层面网格生成的层序分析方法，其特征在于，所述步骤B中，定义一种序，认为两个面片间的序值之差与两面片对应的间距等同。根据此定义，步骤B中构造一目标函数，其含义为：理想情况下，即同相轴追踪结果精确，所有面片的序值正确，同属同一反射界面的面片的序值相等，则序号差值与平均层厚度相等，故目标函数值等于0。一般情况下目标函数值刻画了当前序值与理想序值的差异的大小。上述所述面向三维层面网格生成的层序分析方法目的在于求使目标函数取到最小值的解，即使目标函数描述的差异最小化。

上述所述面向三维层面网格生成的层序分析方法，其特征在于，所述步骤C中，先对部分未知数预赋值。预赋值是指：根据面片的覆盖关系，将代表最上层面片的未知数预先给定一个序值，以约束方程，使得方程有解。该预赋值可以是对序值的某种预期，也可以是相对于某一个参考层面模型的平均距离。然后根据多元函数求极值的方法，将目标函数对每一个没有预赋值的未知数求偏导数，并令结果等于0，处理的结果为得到一线性方程组。该方程组的解即为使目标函数取最小值的解。求解使用Gauss-Seidel迭代法，其收敛性已被证明是有保证的，但不限于Gauss-Seidel迭代法，证明的结果表明凡收敛充分条件为系数矩阵为正定矩阵的求解方法均可。

上述所述面向三维层面网格生成的层序分析方法，其特征在于，所述步骤D中，根据求出的序值给面片或线段分类。该序值具有如下特点：相互分离但是同属一层的面片，其序值接近。从而可以根据该特征将分布在三维空间相互分离但同属于同一层的面片进行合并。为此有两种处理方法，一种是在竖直方向上直接采样；另一种是将序值作为面片的属性在竖直方向上插值形成一个属性体，然后在属性体中追踪指定序值的等值面，然后取竖直方向上分布在等值面附近的面片。分类归并后的结果作为地层建模的输入数据。

附图说明

以下结合附图及具体实施例对本发明进行详细描述，但不作为对本发明的限定。

图1为一测线地震剖面；

图2为图1所示地震剖面的同相轴自动追踪结果，可看成是三维同相轴追踪所得面片集合与该测线剖面的交线；

图3为竖直扫描线对剖面上同相轴进行扫描的示意图；

图4为扫描所得序关系表示为有向图；

图5为同相轴发生窜层的示意图；

图6为同相轴自动追踪线集的层序分析结果图，其中序值映射为灰度值，可以通过灰色与白色过渡产生的分界线看出排序后线段之间的大致关系，反映出了地层的走势；

图7为三维的面片的层序分析结果，其中序值映射为灰度值，通过灰度渐变的效果可以清楚的看出地层面的三维分布；

图8为抽取了某一序值附近的面片的离散点生成的层面网格，黑色为抽取的离散点；

图9是抽取了多个序值生成的多层层面网格；

图10为面向三维层面网格生成的层序分析方法流程图；

具体实施方式

本实例详细介绍一种面向三维层面网格生成的层序分析方法。本发明中，采用和传统同相轴自动追踪不同的思想，转而利用同相轴自动追踪结果携带的信息，从全局出发，通过对面片集合进行层序分析，恢复追踪结果的空间拓扑信息，并利用该信息对结果进行分类，从而用于地质建模。

在本发明中，输入数据为同相轴自动追踪生成的结果集。同相轴的自动追踪可以看成是给定一些种子点，根据种子点的特征在三维空间中搜索临近的性质相似的采样点，不断扩大连通域的过程，其结果为由采样点构成的面片(三维)或线段(二维)。本发明不限定同相轴自动追踪使用的算法，具体可以参阅相关文献。

假设已通过同相轴自动追踪获得一初始的面片集合(步骤1000)，则该面片集合与某一测线剖面相交构成一线集，若整个地震数据体有k条测线，则面片集合可以看作是由k个测线剖面上的线集构成的，而每一个面片则是由分散在多个测线剖面上的线段组成的。图1，图2展示了一测线地震剖面及其自动追踪结果样例，是面片集合与某一测线剖面的交线。故可以以二维剖面为单位处理，逐测线剖面进行。这里将图1，图2抽象为图3的形式，并以图3为例详细说明在二维剖面上实施的具体细节。

本发明的核心思路为：已知连续的同相轴代表了地下的反射界面，即地层分界面。地层在形成的过程中，新的地层总是覆盖在老的地层之上，即使后来发生了断裂，挤压，扭曲等作用，这种已经形成的序关系是稳定的，可以通过这种序关系来确定面片与面片之间的关系。面片层序分析就是给这些面片(线段)排序，并且希望排好序后的序值有如下特点：1)序值的大小符合面片的空间关系，即一个面片的序值始终比被压在它下面的面片的序值小(大)；2)如果两个分离的面片实际上代表的是同一个地层面，那么在排序后，希望它们的序号值之差尽可能的小；3)允许相等的序号。如果有满足以上特点的排序，就可以通过一面片的序号来推断与其它面片的关系。以下都认为上层的序号比覆盖在其下的层的序号小，但不限于此顺序。

根据上述思路，第一步，需要建立包含面片空间拓扑信息的关系式。这里采用的方法是：在竖直方向上对面片进行扫描(步骤1030)。本实施中，扫描线的水平位置是按照交叉线(cross line)的位置取的，考虑到运算量，可以选择降低扫描的密度。这种竖直方向上的扫描可以确定线段(面片)间在竖直方向上的关系。以图3为例，通过竖直方向的扫面线1，可以确定线段A覆盖在线段D之上，根据上述提到的序值特点，设x_A、x_D为线段A、D的序值，则有x_A＜x_D。同样对于扫描线2，可以得到x_A＜x_B＜x_C＜x_D＜x_G，x_B＜x_C＜x_D＜x_G，x_C＜x_D＜x_G，x_D＜x_G。其它扫描线类似。这些扫描得到关系可以表示成有向图的形式，如图4所示，其中箭头指向覆盖的线段。这种关系定性的反映了线段(面片)之间的覆盖关系，还需要加入可以定量的信息。为此引入平均层厚度和序关系计数的概念。

序关系计数是指某一序关系在层序扫描后出现的次数。以图3为例，考察线段A覆盖在别的线段之上的次数，对于全部四条扫描线，则x_A＜x_B出现2次，x_A＜x_C出现1次，x_A＜x_D出现2次，x_A＜x_E出现1次，x_A＜x_G出现2次，对应的序关系计数分别为2、1、2、1、2。这里把序关系计数保存成矩阵的形式，设序关系计数矩阵为R，则对于图3所示的例子，有：

$R=\left(\begin{array}{ccccccc}0& 2& 1& 2& 1& 0& 2\\ 0& 0& 1& 1& 2& 1& 2\\ 0& 0& 0& 1& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right)$

其中矩阵元素R(i，j)表示扫描过程中关系x_i＜x_j出现的次数。这里i≠j，即不考虑对角元素。

在扫描的过程中，除了可以发现线段之间的覆盖关系，还可以得到在该水平位置处，两线段在竖直方向上的间距，如图3所示，扫描线与各线段存在交点，利用这些交点计算两线段之间的间距，并将每一条扫描线得到的间距相加，记录在累积间距矩阵D中。矩阵D的维数和R相同，矩阵元素D(i，j)表示对应于关系x_i＜x_j，线段i与线段j在扫描过程中得到的间距的累加和。和R一样，D的对角元素也为0，后续运算中不使用对角元素。由此可知这种间距是与序相关的，D(i，j)与D(j，i)的含义并不相同，其累加和表示扫描中必须出现对应的关系。利用累积间距矩阵D，可以得到平均层厚度的概念。设平均层厚度则有：

$\overline{d_{\mathrm{ij}}}=\frac{D(i,j)}{R(i,j)}$

易知平均层厚度也是序相关的，这里i≠j。

利用平均层厚度和序关系计数，我们使用最小二乘法，从全局的角度，试图给每一个面片赋予一个序值，该序值满足我们的某个优化目标。这里我们给出的目标函数如下：

$F\left(\overrightarrow{x}\right)=\underset{i=1}{\overset{n}{\Sigma }}\underset{\underset{j\ne i}{j=1}}{\overset{n}{\Sigma }}{(x_j-x_i-{\bar{d}}_{\mathrm{ij}})}^2·C_{\mathrm{ij}}^2$

该目标函数的含义是：把两面片的序号差值x_j-x_i与平均层厚度等同起来，认为理想情况下，即同相轴追踪结果精确，所有面片的序值正确，同属同一反射界面的面片的序值相等，则序号差值x_j-x_i与平均层厚度相等，故目标函数等于0。从而我们的目标就是求解目标函数的一解向量该解向量使得的值最小。解向量的各分量值即为对应面片的序值。目标函数中的作用是控制解的倾向性。是序关系的可信度，越大，说明在扫描的过程中关系x_i＜x_j出现的越多，若x_i和x_j的误差越大，则对的影响越大，起到了惩罚的作用。同时也能压制自动追踪错误(如窜层，见图5)带来的影响。

通过对数据体逐剖面的扫描，生成矩阵R和D，可以构建出我们需要的目标函数。第二步要求出使目标函数达到最小值的解。目标函数无法直接求解，需要预处理。首先要对一部分的未知数赋初值，原因在于目标函数中只有计算相对差值的项，无法将解固定，需要消去一部分未知数，方法就是使目标函数中的部分x_i变成常数值。对于预赋值x_i的选择，这里采用了图论的相关知识。已知扫描的结果可以表示为有向图的形式，如图4所示，节点A是图中入度为0的节点，通过对节点A预赋值，则可以通过传导约束其余的节点。节点A的含义是面片中处于最上层的面片，没有其它面片覆盖在其之上，在关系矩阵R中表现为该面片对应的列向量为0向量。根据此特点，较容易从未知数中选出需要预赋值的未知数。由于这些面片处于顶部，预赋值的内容可以选择为这些面片到顶边界的平均距离。也可以根据地层的起伏趋势，在这些面片上方外加一个吻合起伏趋势的参考层面，预赋值为面片到参考层面的平均距离，效果更好。

对目标函数预赋值后，根据多元函数求极值的方法，将目标函数对每一个未赋值的未知数求偏导，可以得到一线性方程组(步骤1040)：

其中，Ψ为没有预赋值的未知数的编号集合，Φ为预赋值的未知数的编号集合，Ψ(i)表示对Ψ的元素按序号升序排序后取第i个元素。方程组的系数矩阵为一稀疏矩阵，适合使用迭代法求解。这里选择使用Gauss-Seidel迭代方法(步骤1050)。在实际使用中，可不限于Gauss-Seidel迭代方法，我们已经证明，经过预赋值处理的目标函数生成的线性方程组的系数矩阵是对称正定阵，凡是以对称正定阵为收敛的充分条件的方法都可以。

第三步，将求得的解作为序值赋给对应的面片。根据序值的特点，将具有较接近序值的面片合并分组，然后作为地质建模的数据(步骤1060)。这里可以使用两种方式，一种是在竖直方向上直接采样。例如有大量面片的序值在x附近，则以x为标准，在竖直方向上采样时(类似前述的扫描过程)，选择具有与x最接近的序值的面片，并且竖直方向上这样的面片只取一个，避免重叠的现象；另一种是将序值作为面片的属性在竖直方向上插值形成一个属性体，然后在属性体中追踪指定序值的等值面，然后取竖直方向上分布在等值面附近的面片。如此完成了同相轴自动追踪结果的分类，将零碎的面片合并成了较完整的层面点集，该点集可直接用于层面网格建模。

另外，根据目标函数中含有的的作用易知排序的结果易受扫描统计出的关系计数影响，由此可以得到对排序的结果施加约束的方法(步骤1010，1020)。由于断层的存在，同相轴自动追踪往往会发生中断，这使得序关系也发生中断，对应于有向图，其特征为有许多分离的子图。为提高排序的精确度，可以在扫描前加入约束层，该约束层应当能够覆盖希望影响的区域，以使得扫描的结果中较大。约束层应当反映正确的地层结构，其来源可以是以人工的方式对地震数据体的有选择的预解释。另外也可以通过人工合并某些面片扩大面片的影响范围，同样能够改变排序的结果。

根据以上步骤，给出如下实施例，如图6到图9所示。

图6为图2剖面同相轴自动追踪线集的分析结果。将序值映射为灰度值，可以通过灰色与白色过渡产生的分界线看出排序后线段之间的大致关系，反映出了地层的走势。对比图1和图2，可见分析的结果与实际的地层空间特征相吻合(观察图1和图2中地层界面的分布趋势)。

使用和图6相同的二维处理方式，通过采用逐剖面处理的方式，扫描完全部的测线剖面后可以得到三维的面片分析结果。对结果进行灰度映射，效果如图7所示。通过灰度渐变的效果可以清楚的看出地层面的三维分布。

图8为抽取了某一序值附近的面片的离散点生成的层面网格，黑色为抽取的离散点。图9是抽取了多个序值生成的多层层面网格。通过实施例可以看出方法的有效性。

本发明还可以有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。