CN104252549A - 一种基于克里金插值的分析布井方法 - Google Patents

Abstract

一种基于克里金插值的分析布井方法，属于计算机分析技术领域。其中包括：分析原始采样点，区域内每口井包含4个属性值（井号，横坐标，纵坐标，高程值）；确定布井区域，默认的布井区域为包含原始数据点的最大矩形；并且在整个布井区域内建立由行列索引的网格；将每个网格以其中一条对角线划分为两个三角形，将整个布井区域建立一个互不交叉的三角形列表；采用克里金插值方法对每个所述网格点插值；等值线追踪，遍历所有三角形，在三角形内寻找各个通过的等值线；填充区域的搜索，在所有的三角形内部，按照通过的等值线将三角形划分为多个区域，确定布井的点。

Description

一种基于克里金插值的分析布井方法

技术领域

本发明涉及一种基于克里金插值的分析布井方法，属于计算机分析技术领域。

背景技术

通常通过地震勘测等先进科学方法能够采集到各种地质数据。但是由于勘测成本极高，地质数据的分布往往是不规则的且离散的，而且在许多情况下数量有限。那么怎样才能通过有限的地质数据来构造出完整的地质模型，并计算出未知区域呢？

早期的地质人员在布井区域中手工绘制采集的数据点，并利用三角网线性插值的方法描绘出所有的等值线。绘制一张标准的等值线图会耗费大量的时间，只能用图纸的方式保存下来。后面采集数据点增加后也无法支持修改，必须重新绘制一张等值图，大量的重复工作需要进行。

现在，通过克里金插值算法，使得采集的原始数据包含的地理特征无明显损失地传递到被估算的网格点，通过有限的地质数据来构造出完整的地质模型。由此获得的油藏地质模型描述更加逼近真实情况，将地下油藏的内部状态形象直观地呈现给油田管理技术人员，使他们做出正确的决策，提高布井的准确性，显著提高工作效率。

克里金插值方法不仅考虑已知点和被估计点的相对位置，还考虑了已知点之间的相对位置关系。克里金分两步，首先对空间场进行结构分析，提出变差函数模型，最后根据变量图函数模型进行克里金计算。克里金插值法是一种线性、无偏、方差最小的空间插值方法。计算时以空间结构分析为基础进行估值，插值过程中可以反映空间场的各向异性，并且充分利用数据点之间的空间相关性。优点：插值数据的真实性；缺点：计算时常常以耗费时间和资源为代价。

发明内容

本发明的目的在于提供一种基于克里金插值的分析布井方法，解决上述问题。

为了实现上述目的，本发明采用的技术方案如下：

一种基于克里金插值的分析布井方法，包括以下步骤：

步骤1，分析原始采样点，区域内每口井包含4个属性值（井号，横坐标，纵坐标，高程值）；

确定布井区域，默认的布井区域为包含原始数据点的最大矩形；并且在整个布井区域内建立由行列索引的网格；

步骤2，在每个网格点上，建立普通克里金方程组，求解网格点的值；主要包括：变差函数模型的选择，搜索策略的选择，普通克里金方程组的建立，方程组的求解，以及根据系数求解网格点的值；

步骤3，将每个网格以其中一条对角线划分为两个三角形，将整个布井区域建立一个互不交叉的三角形列表；

步骤4，基于三角形的追踪原则，与矩形网格追踪相比较，在三角形的三边上只可能有两个相同的等值点，这样在一个三角形内只可能有唯一的一条等值线，编程实现简单；而在矩形内有可能有四个相同的等值点，这样就有内部可以形成多种组合形式，造成编程实现复杂；

步骤5，在所有的三角形内部，按照通过的等值线将三角形划分为多个区域；

步骤6，单个三角形内部区域搜索：单个三角形内部通过的等值线只有两种情况，要么没有任何等值线通过，要么有一条或多条等值线通过；

为确定被等值线分割的三角形区域，每个区域只需要找到对应的边界点即可；和区域相关的点的集合仅仅由三角形的三个顶点以及等值线与三角形边的交点组成；将点集合中的每个点都设置两个关键字，即每个点都包含双关键字：高程值和顺序标志；其中顺序标志默认为空；确定布井的点。

附图说明

当结合附图考虑时，通过参照下面的详细描述，能够更完整更好地理解本发明以及容易得知其中许多伴随的优点，但此处所说明的附图用来提供对本发明的进一步理解，构成本发明的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定，如图其中：

图1原始采集点的分布图；

图2网格化和三角化后形成的网格点分布；

图3单个三角形内部通过的等值线情况示意图；

图4等值线追踪过程中非封闭等值线追踪过程示意图；

图5等值线追踪过程中封闭等值线追踪过程示意图；

图6单个三角形内部搜索等值区域示意图；

图7等值线图最终效果图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

如图1、图2、图3、图4、图5、图6、图7所示，一种基于克里金插值的分析布井方法，含有以下步骤；

1.克里金插值算法：

分析原始采样点：如图1，区块内每口井包含4个属性值（井号，横坐标，纵坐标，高程值）；

网格线索几何化：获取所有井坐标中的最大值，最小值，计算一个矩形区域用以确定区块范围。根据区块矩形计算规则网格点的坐标值；

变差函数模型的选择：变差函数模型是指以空间两点间距离为自变量，具有解析表达式的函数。变差函数是克里金插值中不可缺少的组成部分，将直接参与克里金插值的估算。常用的变差函数模型：线性模型，球状模型等等。线性模型：γ(h)=C×h;其中h表示两点的距离，C表示比例值（在应用软件中，C默认等于1，同时支持用户对上述参数的自定义设置），γ(h)函数值会参与之后的克里金估算。球状模型： $\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ C_0+C(\frac{3}{2}\&times;\frac{h}{a}-\frac{1}{2}\&times;\frac{h^3}{a^3})& 0<h\&le;a\\ C_0+C& h>a\end{array}\right.;$ 其中h表示两点的距离，C₀表示块金常数，C表示拱高，a表示变程。（在应用软件中默认C₀=1，C=1，a=10，同时支持用户对上述参数的自定义设置）；

搜索策略:为了计算网格点位置x₀的值Z(x₀)，那么选取哪些原始采样点来计算它呢？通过搜索策略来寻找，假设搜索出n个与待估点相关的采样点,设为Z(x_i),i=1,2...n。利用Z(x_i)的线性组合来估算Z(x₀)。

通常情况下，都采用全部搜索，即选取所有的采样点加入计算，这样会保证计算的真实性。但是如果原始采样点数量非常大，如果利用全部搜索，计算量会直线上升，计算花费的时间将成为考虑的主要因素。那么采用其它搜索策略（如四分搜索）来计算，以提高效率，当然此时的插值的真实性将会受到影响。质量与效率往往不可兼得，用户通过自定义的操作来选择采用更折中的方案;

普通克里金方程： $\left\{\begin{array}{cc}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{i}C(x_i,x_j)-\mathrm{\&mu;}=C(x_0,x_j)& (j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n)\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_i=1& \end{array}\right.;$ $\left\{\begin{array}{cc}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{i}C(x_i,x_j)-\mathrm{\&mu;}=C(x_0,x_j)& (j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n)\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_i=1& \end{array}\right.;$

其中C(x_i,x_j)=γ(x_i,x_j)表示变差函数计算来的函数值，μ表示拉格朗日（Lagrange）乘子，λ_i表示加权系数，x₀表示被估算点，n表示参与计算的采样点。

以上（n+1）个方程的方程组对应的矩阵形式： $\left[\begin{array}{ccccc}{C}_{11}& C_{12}& ...& C_{1n}& 1\\ C_{21}& C_{22}& ...& C_{2n}& 1\\ .& .& & .& .\\ .& .& ...& .& .\\ .& .& & .& .\\ C_{n1}& C_{n2}& ...& C_{\mathrm{nn}}& 1\\ 1& 1& ...& 1& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\&lambda;}}_1\\ {\mathrm{\&lambda;}}_2\\ .\\ .\\ .\\ {\mathrm{\&lambda;}}_n\\ -\mathrm{\&mu;}\end{array}\right]=\left[\begin{array}{c}{C}_{01}\\ C_{02}\\ .\\ .\\ .\\ C_{0n}\\ 1\end{array}\right];$ 其中C_ij=C(x_i,x_j)；

方程求解：利用高斯消元法求解方程（LU矩阵分解），求解向量 $\stackrel{\&RightArrow;}{\mathrm{\&lambda;}}=\left[\begin{array}{cccc}{\mathrm{\&lambda;}}_1& {\mathrm{\&lambda;}}_2& \&CenterDot;\&CenterDot;\&CenterDot;& {\mathrm{\&lambda;}}_n\end{array}\right];$

计算被估算点：通过与Z(x_i)的线性组合来计算被估算点， $x_0=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_i\&CenterDot;Z\left(x_i\right)={\mathrm{\&lambda;}}_1\&CenterDot;Z\left(x_1\right)+{\mathrm{\&lambda;}}_2\&CenterDot;Z\left(x_2\right)+\&CenterDot;\&CenterDot;\&CenterDot;+{\mathrm{\&lambda;}}_n\&CenterDot;Z\left(x_n\right);$ 图2为最后计算出的规则网格格点下的值。

2.三角化：

三角化的目标是把网格化处理后生成一系列的三角形列表，为等值线的追踪算法和等值线的填充区域搜索算法做准备。

通过网格化，离散的测量采集数据点可以转化为连续的网格数据曲面。由于网格化是规则的矩形行列分布，所以对等值线的追踪算法以及等值线的填充区域搜索算法自然地联想到基于四边形来完成。

但是考虑到四边形追踪涉及的情况复杂多变，不利于解决问题。

三角形相对于四边形来说，顶点和边都要更少，基于三角形的追踪更加简单，复杂问题简单化了。

网格化后的离散点以规则的矩形分布，每个矩形内部可以划分为右上和左下两个三角形，这样可以确定三角化的唯一性。

生成的三角性列表可以方便地从左到右或从上到下访问到，也满足了格网数据的规则性。最后三角化处理后的数据如图2。

三角化最终的数据处理结果为：网格点列表，边列表，三角形列表。点(网格点)：包含实际坐标值(x,y)；点对应的高程值(elevation)等等。边：一条边包含2个点，并且最多所属2个相关三角形。

数据结构如下：边的类型(边界边，非边界边)；2个端点(点列表所属的索引值p1,p2);所属的2个三角形(triangle1,triangle2)。

三角形：一个三角形包含3个点，同时包含了3条边。

数据结构如下：3个顶点(p1,p2,p3)。点，边，三角形内部都包含了丰富的几何要素，对程序内部处理时，都能够方便地获取到它关联的几何要素。

等值线追踪算法：

等值线追踪的目标是以确定的高程值为依据，搜索该等值线与三角形边相交的所有交点，最后形成等值线相关的控制点列表。

等值线的追踪算法完全基于三角形，在三角形列表中采用完全搜索。

(1)预准备：网格点列表；边列表；三角形列表。如果边与等值线相交，那么称此边为等值边。

(2)非闭合等值线的搜索：

①根据等值线的高程值，在边界边中搜索一条等值边作为起始边，如果不存在，则说明该等值线闭合，执行(3)。并设置该边为已搜索标志。通过线性插值计算起始边与等值线的交点坐标值，并记录为该等值线的控制点。

②在起始边相关三角形的其它边上寻找第二条等值边，找到后线性插值出控制点坐标，并设置该边为已搜索标志。

③此三角形设置为已搜索标志，则当前的等值边只有一个相关三角形，在该三角形内搜索下一条等值边。

④递归执行③，直到找到的最后一条等值边是边界边。

(3)闭合等值线的搜索：

①若在所有的非边界边中都无法找到起始边，则说明该等值线无法分析，应退出该等值线的分析而开始另一条等值线的搜索。

②在起始边的一个相关三角形中寻找第二条等值边，找到后线性插值出控制点坐标，并设置该边位已搜索标志。

③等值线已经过的三角形设置为已搜索标志，则当前的等值边只有一个相关三角形，在该三角形内搜索下一条等值边。

④递归执行③，直到找到最后一条等值边是起始边。

(4)对于所有高程值，循环执行(2)、(3)直到找到所有等值线。

基于三角形的等值区域搜索；

等值线的填充区域搜索可以使区域显示呈现阶梯式的颜色渐进变换，同时也可以突出的显示重点区域。本系统提出了一种基于三角形的算法。

预准备：已知所有的高程值数组；已知三角形列表；

单个三角形内部区域搜索：单个三角形内部通过的等值线只有两种情况，要么没有任何等值线通过，要么有一条或多条等值线通过。假设该三角形的三个顶点对应的高程值分别为20（小）,50（中）,100（大）。

遍历高程值数组，如果高程值位于（20,50）或（50,100），则说明该高程值肯定通过此三角形。否则没有任何等值线通过。

假设如图6中，高程值45的等值线通过（20,50）的三角形边，高程值65和85通过（50,100）的三角形边。为确定被等值线分割的三角形区域，每个区域只需要找到对应的边界点即可。和区域相关的点的集合仅仅由三角形的三个顶点以及等值线与三角形边的交点组成。将点集合中的每个点都设置两个关键字，即每个点都包含双关键字：高程值和顺序标志。其中顺序标志默认为空。

如图7中，15个点包含的关键字经过双关键字排序后如下：（20，45(1)，45(2)，45(3)，45(4)，50,65(1)，65(2)，65(3)，65(4)，85(1)，85(2)，85(3)，85(4)，100）。因为通过的等值线的条数为n，那么划分的区域为n+1。故将上述的点列表进行区域搜索：20，45(1)，45(2)||45(3)，45(4)，50,65(1)，65(2)||65(3)，65(4)，85(1)，85(2)||85(3)，85(4)，100||。

如果三角形内部没有等值线通过，那么有三角形三个顶点围成的区域记为等值线区域。

循环三角形列表中的所有三角形，重复步骤2即可完成对所有区域的搜索。

如上所述，对本发明的实施例进行了详细地说明，但是只要实质上没有脱离本发明的发明点及效果可以有很多的变形，这对本领域的技术人员来说是显而易见的。因此，这样的变形例也全部包含在本发明的保护范围之内。

Claims (1)

1.一种基于克里金插值的分析布井方法，其特征在于包括如下步骤：

步骤1，分析原始采样点，区域内每口井包含4个属性值（井号，横坐标，纵坐标，高程值）；

确定布井区域，默认的布井区域为包含原始数据点的最大矩形；并且在整个布井区域内建立由行列索引的网格；

步骤2，在每个网格点上，建立普通克里金方程组，求解网格点的值；主要包括：变差函数模型的选择，搜索策略的选择，普通克里金方程组的建立，方程组的求解，以及根据系数求解网格点的值；

步骤3，将每个网格以其中一条对角线划分为两个三角形，将整个布井区域建立一个互不交叉的三角形列表；

步骤4，基于三角形的追踪原则，与矩形网格追踪相比较，在三角形的三边上只可能有两个相同的等值点，这样在一个三角形内只可能有唯一的一条等值线，编程实现简单；而在矩形内有可能有四个相同的等值点，这样就有内部可以形成多种组合形式，造成编程实现复杂；

步骤5，在所有的三角形内部，按照通过的等值线将三角形划分为多个区域；

步骤6，单个三角形内部区域搜索：单个三角形内部通过的等值线只有两种情况，要么没有任何等值线通过，要么有一条或多条等值线通过；

为确定被等值线分割的三角形区域，每个区域只需要找到对应的边界点即可；和区域相关的点的集合仅仅由三角形的三个顶点以及等值线与三角形边的交点组成；将点集合中的每个点都设置两个关键字，即每个点都包含双关键字：高程值和顺序标志；其中顺序标志默认为空；确定布井的点。