CN102982566B - 一种基于最小曲率法插值的含断层的等值线图自动生成方法 - Google Patents

Abstract

一种基于最小曲率法插值的含断层的等值线图自动生成方法，包括以下步骤：1，网格化，确定研究区域，并在所述研究区域内建立网格；2，三角化，将每个所述网格以其中一条对角线划分为两个三角形，并将整个所述研究区域建立一个互不交叉的三角形列表；3，数据的预处理，建立点列表、边列表和三角形列表，同时记录边类型；4，采用最小曲率法插值方法对每个所述网格点插值；5，采用虚拟点技术进行所述等值线追踪；6，填充区域的搜索；7，等值线数据的格式化，所述等值线数据生成后，将所述等值线数据和所述填充区域的数据格式化后保存。本发明的等值线图的生成方法，使等值线图生成速度快，停靠精确，出图效果美观，有利于对整个等值区域的研究。

Description

一种基于最小曲率法插值的含断层的等值线图自动生成方法

技术领域

本发明涉及等值线图自动生成技术领域，尤其涉及一种基于最小曲率法插值的含断层的等值线图自动生成方法。

背景技术

等值线图是一类在地质勘探、油气田开发、气象研究等领域应用极其广泛的重要图件。以油气田地质勘探为例，通常通过地震勘测等先进科学方法能够采集到各种地质数据，但是由于勘测成本极高、数据采集困难等因素，导致采集的地质数据数量有限，而且分布往往是离散且不规则的，不利于分析整个区域。但是，通过插值技术，能够根据有限的地质数据，计算出整个区域内所有点的数据，从而将采集的离散数据转化为连续的曲面数据，并构建出完整的地质模型，将地下油藏的内部状态形象直观地展现出来。

目前业界生成等值线图的基本流程分为以下几个步骤：

(1)网格化：采用矩形网格化，建立行列索引的网格覆盖整个研究区域。

(2)网格点的插值方法：采用距离反比加权插值或克里金插值等方法。

(3)等值线的追踪方法：采用基于四边形的追踪方法。

(4)填充区域的搜索方法：采用三角形内部的等值区域填充方法。

目前业界实现等值线的技术方法如下：

(1)如何插值：应用最为广泛的是克里金插值方法。它是一种线性、无偏、方差最小的空间插值方法，以空间结构分析为基础进行估值，充分利用了数据空间场的性质，在插值过程中可以反映空间场的各向异性，并且充分利用数据点之间的空间相关性，是一种最优内插法。

(2)如何追踪：通过基于四边形的追踪方法对等值线进行追踪。在网格化后的每一个矩形内线性插值计算等值点，并以此等值点为起点在相邻的矩形中追踪剩余的等值点。

(3)如何停靠到断层线：目前国内外采用的断层停靠方法多样，如：

白化技术：等值线不停靠到断层线，而是沿着断层线边缘进行追踪；

垂直停靠技术：当等值线到达断层边缘附近时，获取等值点在断层线上的投影点作为停靠点。

(4)如何填充：通过遍历每个三角形内部的等值区域进行填充，达到渐变颜色的无缝过渡。

以上实现方法各自的缺点及原因：

(1)首先，克里金插值方法在计算每一个网格点时部要重新计算原始数据点的权值，即每次计算一个网格点都要解一个n元一次方程组，故耗费的时间以及内存部是非常大的，因此克里金插值方法效率上有很大的局限性。其次，克里金插值方法是在整个区域连续的假设下提出的，但是加入断层后，断层破坏了区域原有的连续性，因此克里金插值方法无法直接应用于带断层的区域。最后，通过建立断层多边形，在每个断层多边形局部使用克里金插值。由于断层分为正、逆断层，且分布也毫无规则，因此准确合理地划分断层多边形非常困难，需要对多种形态的断层分别处理，增加了问题的复杂程度。

(2)基于四边形的等值线追踪方法。在处理一个网格内同时出现四个等值点的情况时，如何选择下一个等值点会很困难。针对这样的情况，需要建立双元双三次一阶导数连续的曲面函数预测下一个等值点的选择。

(3)如果采用白化技术，等值线没有停靠到断层线上，而是沿着断层的边缘追踪，这样会出现断层线边缘附近会有多条等值线聚集，仅仅采用屏幕白化遮挡了断层线边缘的等值线，达到等值线停靠断层的假象，导致等值线图效果往往不美观；如果采用垂直停靠技术，以等值点在断层线上的投影作为停靠点，这样会出现等值线停靠到断层外、等值线出现弯折等停靠不精确、不美观的情况。

(4)每个三角形内部可能有多条等值线通过，这些等值线把三角形划分出了多个小等值区域，再根据等值线的高程值确定这些小等值区域的值，通过不同的颜色进行填充，但由于小的等值区域在屏幕显示时非常小，无法进行纹理填充。

发明内容

本发明的目的在于设计一种新型的基于最小曲率法插值的含断层的等值线图自动生成方法，解决上述问题。

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

一种基于最小曲率法插值的含断层的等值线图自动生成方法，包括以下步骤：

步骤1，网格化，确定研究区域，并在所述研究区域内建立网格；

步骤2，三角化，将每个所述网格以其中一条对角线划分为两个三角形，并将整个所述研究区域建立一个互不交叉的三角形列表；

步骤3，数据的预处理，建立点列表、边列表和三角形列表，同时记录边类型；

步骤4，网格点的插值，采用最小曲率法插值方法对每个所述网格点插值；

步骤5，等值线追踪，采用虚拟点技术进行所述等值线追踪；

步骤6，填充区域的搜索，建立基于边框搜索线和断层搜索线形成的环结构，按照逆时针方向遍历搜索线寻找所述填充区域；

步骤7，等值线数据的格式化，所述等值线数据生成后，将所述等值线数据和所述填充区域的数据格式化后保存。

优选的，所述研究区域为包含原始数据点和断层线的最大矩形，所述原始数据点包含的属性有横坐标x、纵坐标y、高程值z。

优选的，所述的边类型包括边框边、内部边、断层边和虚拟边，所述边框边是指边位于边框上；所述断层边是指边有断层通过；所述虚拟边是指边本身不存在，随着等值线的追踪动态创建；所述内部边是指除所述边框边、所述断层边和所述虚拟边的边。

优选的，所述网格点插值主要包括：

(1)所述网格点的迭代方程的建立；

(2)所述网格点的初始值设置；

(3)所述网格点的迭代方程的迭代求解。

优选的，建立所述网格点的迭代方程时，首先建立以所述网格点为中心的中心区域，在所述中心区域内，通过将多个所述原始数据点加权平均，使得每个所述网格点的矩形区域至多有一个原始数据点与之对应。

优选的，不同所述网格点的迭代方程包括：分别是远离边框的网格点迭代方程、边框外的虚拟点的迭代方程、靠近断层的虚拟点的迭代方程、受原始数据点约束的网格点的迭代方程和断层端点附近的网格点的迭代方程。

优选的，在建立所述断层端点附近的网格点的迭代方程时，通过周围所述网格点加权平均来设置。

优选的，建立边框附近的所述网格点的迭代方程和断层附近的所述网格点迭代方程时，当出现邻近网格点不足的情况时，附加两排虚拟网格点，通过所述虚拟网格点所属的所述网格点的迭代加权平均来建立所述虚拟网格点的迭代方程。

优选的，所述网格点的初始值的设置方法如下：以受所述原始数据点约束的所述网格点为种子，向上、下、左、右以及左上、右上、左下、右下八个方向扩散，在扩散过程中查看对应的边是否为断层边，如果为断层边，则说明被断层阻断，朝此方向的扩散结束；循环每个种子点的扩散，每个所述网格点都会得到一个关于所述种子点和扩散次数的列表，根据所述列表，通过反距离加权法计算所述网格点的初始值。

优选的，所述等值线追踪采用基于三角形的基本原理，同时采用所述虚拟网格点技术，拓展断层单侧数据的连续性，在所述等值线追踪过程中，时刻检查所述等值线是否和断层相交，并选择所述等值线和断层相交的点做为停靠点。

本发明的有益效果可以总结如下：

1.等值线图能够快速生成；

2.等值线追踪算法更加简单、容易实现；

3.等值线的追踪算法基于三角形，在三角形中追踪等值点可以避免在四边形追踪中会出现选择困难的问题；

4.等值线能够精确地停靠到断层线；

5.等值区域的纹理填充基于边框或断层的搜索线，能够找到完整的等值区域，避免了其它方法中无法进行纹理填充的问题。

附图说明

图1是本发明的等值线图自动生成技术的流程图；

图2是本发明的网格化示意图；

图3是本发明的三角化示意图；

图4是本发明以网格点为中心的区域示意图；

图5是本发明远离边框的网格点与邻近点的位置关系示意图；

图6是本发明靠近边框的网格点与邻近点的位置关系示意图；

图7是本发明黑色矩形所处的网格点在断层的附近的示意图；

图8是本发明黑色矩形所处的网格点为断层附近的虚拟点的示意图；

图9是本发明受原始数据点约束的网格点示意图；

图10是本发明断层端点附近的网格点示意图；

图11是本发明断层端点附近的网格点示意图；

图12是本发明等值线循环追踪示意图；

图13是本发明等值线停靠到边框线的追踪示意图；

图14是本发明等值线停靠到断层线的追踪示意图；

图15是本发明基于边框搜索线和断层搜索线形成的环结构示意图；

图16是本发明由封闭的等值线形成的填充区域示意图；

图17是本发明由非封闭的等形成的填充区域示意图；

图18是本发明等值线图最终生成的总体效果图；

图19是本发明等值线的断层停靠效果图；

图20是本发明等值区域的纹理填充效果图。

具体实施方式

为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例一：

如图1所示的本发明的等值线图自动生成技术的流程图，包括以下步骤：

1.网格化：确定研究区域，默认的研究区域为包含原始数据点、断层线的最大矩形。并且在整个研究区域内建立由行列索引的网格。

2.三角化：将每个网格以其中一条对角线划分为两个三角形，将整个研究区域建立一个互不交叉的三角形列表。

3.数据的预处理：建立点列表，边列表，三角形列表。同时记录边类型：边框边、内部边、断层边、虚拟边。

4.网格点的插值：采用最小曲率法插值方法对每个网格点插值。其主要内容包括：网格点的迭代方程的建立、网格点的初始值设置、迭代求解。

5.等值线追踪：基于三角形的追踪原则，采用虚拟点技术进行等值线追踪。

6.填充区域的搜索：建立基于边框搜索线和断层搜索线形成的环结构，按照逆时针方向遍历搜索线寻找填充区域。

7.等值线数据的格式化：等值线数据生成后，为提高程序的执行效率，将等值线以及填充区域数据格式化，可以入库或放入文件中保存，供以后使用。

本发明采用最小曲率法插值，并应用了虚拟点技术，虚拟点技术拓展了断层单侧数据的连续性，使得最小曲率法具备直接应用断层线约束的能力。该插值方法不但速度极高，而且随着原始数据点的增加效率不受影响，特别适合于研究大量的原始数据点的区域，与其它插值方法相比，避免了插值不准确、效率低的问题。

本发明采用基于三角形的等值线追踪方法，三角形是最简单的多边形，如果等值线进入其中一条边，那么其它两条边只有一个出口，逻辑处理更加简单，容易实现。与基于四边形追踪的方法相比，避免了在四条边上都有等值点的情况下判断困难的问题。

本发明由于应用虚拟点技术，拓展了断层单侧数据的连续性，等值线会外延到断层线另外一侧虚拟点覆盖的区域。在等值线追踪过程中，会时刻检查等值线是否和断层相交，只有和断层相交时才会选择此交点为停靠点，这样的停靠方法精确、自然、效果美观。与其它停靠方法相比，避免了停靠不精确、出图效果不美观的问题。

本发明采用基于边框和断层的区域搜索方法。边框线采用逆时针方向，而每一条断层线的两侧设定两个相反方向，从而能够将每个等值区域形成一个闭合的环。按照有向线段设定的方向搜索，获得完整的等值区域。这样，在这些等值区域内，可以方便地进行颜色填充或纹理填充。不仅实现了纹理填充的效果，而且更有利于研究整个等值区域。

实施例二：

进一步说明本发明的基于最小曲率法插值的含断层的等值线图自动生成技术的具体方法。

第一步，网格化。如图2所示的本发明的网格化示意图，小圆圈表示原始数据点，通常情况下，这些原始数据点来自地震勘探解释数据，分布不均且数量有限。其中原始数据点包含的属性有横坐标x、纵坐标y、高程值z。最外面的矩形确定的区域为研究区域，研究区域的四条边叫做边框。将整个研究区域划分为若干紧密排列的矩形网格，划分为10行和10列(行数和列数可以由用户自定义)，行索引从左至右依次增大，列索引从下到上依次增大，那么整个研究区域被划分为9×9个网格区域，内有10×10个网格点。利用有限的已知数据，计算其它未知区域，通过仅有的7个原始数据点，插值获取到100个网格点的值，则将离散的测量数据转化为连续的数据曲面，有利于对问题的整体研究。

第二步，三角化。如图3所示的本发明的三角化示意图，将每个矩形网格区域内，连接左上角和右下角的对角线，由此形成了左下角和右上角两个三角形；整个研究区域内形成了多条对角线a。同时在每个矩形网格区域内，连接左下角和右上角的对角线，这样整个研究区域内形成了多条对角线b。

通过对角线a，研究区域内形成一个互不交叉的三角形列表，这些三角形列表主要用于等值线的追踪。对角线b形成了多条边，这些边将用于插值计算，但是不会参与等值线的追踪。

第三步，数据预处理。网格点的插值以及等值线的追踪需要很多的数据支持，需要提供数据包括网格点列表、边列表、三角形列表。

网格点列表：所有网格点组成的集合，每个网格点包含行索引号i，列索引号j，虚拟点列表，在断层附近或边框附近的网格点需要虚拟点，故建立一个虚拟点列表。

边列表：每个矩形网格区域内部有四条边，以及两条对角线，这些所有的边组成的集合叫做边列表。每条边包含边类型，边上的两个端点对应的网格点索引，边所属三角形(最多两个)。边类型分为四种，分别是：边框边，内部边，断层边，虚拟边。其中边位于边框上，则称为边框边；如果边有断层通过，则称为断层边；如果边本身不存在，是随着等值线的追踪动态创建的，则此边称为虚拟边；否则，其它部称为内部边。如果边为断层边，需要记录通过它的断层号索引，断层点的索引。如果边是虚拟边，需要记录它所属的网格点。

三角线列表：每个矩形网格被红色的对角线分为两个三角形，整个研究区域内所有的这些三角形组成的集合叫做三角形列表，每个三角形包含三个顶点，三条边。

第四步，网格点插值，包括以下内容：

建立以网格点为中心的区域：

如图4所示的本发明以网格点为中心的区域，原始数据点不规则地散落在整个研究区域内，每个网格点受所有原始数据点的共同约束。但是为了问题的简化，特别是建立网格点的迭代方程时，一个网格点最多仅受一个原始数据点约束。

虚线矩形框代表了对应的以网格点为中心的矩形区域，而包含在此区域的原始数据点可能有0个或1个甚至多个。为了迭代方程的建立，每个网格点的矩形区域至多有一个原始数据点与之对应。那么多于一个原始数据点，则将多个原始数据点加权平均

$(\frac{x_1+x_2+x_3}{3},\frac{y_1+y_2+y_3}{3},\frac{z_1+z_2+z_3}{3})$

即相当于取所有的原始数据点的重心。

(2)远离边框的网格点迭代方程

如图5所示的本发明远离边框的网格点与邻近点的位置关系示意图，黑色矩形所处的网格点为中心点，圆圈所处的网格点称为内环点，矩形所处的网格点称为中环点，五角星所处的网格点称为外环点。中心点受内环点，中环点，外环点约束，并且约束的条件越来越弱。通过中心有限差分近似得到：

$z_{00}=\frac{\left\{\begin{array}{c}[4(1+{\mathrm{\α}}^2)(1-T_I)+T_I\left]\right[z_{10}+z_{-10}+{\mathrm{\α}}^2(z_{01}+z_{0-1})]\\ -(1-T_I)[z_{20}+z_{-20}+{\mathrm{\α}}^4(z_{02}+z_{0-2})+2{\mathrm{\α}}^2(z_{11}+z_{-11}+z_{1-1}+z_{-1-1})]\end{array}\right\}}{[(6+8{\mathrm{\α}}^2+6{\mathrm{\α}}^4)(1-T_I)+2(1+{\mathrm{\α}}^2)T_I]}$ (公式1)

其中α为网格的横纵比Δx为网格矩形的长，Δy为网格矩形的宽，T_I为内部张力(默认为0)。

(3)边框外的虚拟点的迭代方程

如图6为靠近边框的网格点与邻近点的位置关系示意图，黑色圆圈表示的网格点部是落入研究区域的真实的网格点，其中由虚线多边形包围的网格点处于边框附近。边框附近的网格点通过公式1在使用时会出现邻近网格点不足的情况，采取的方法是在边框外附加两排虚拟的网格点，称为虚拟点。这样所有的网格点部能够根据公式1建立迭代方程。

附加的虚拟点根据边界条件建立差分方程，下面的公式仅仅是x方向上的边框附近的虚拟点的迭代方程；y方向的可以类似获取。

框外的第一层虚拟点的迭代方程：

$z_{-10}=\frac{2(1-T_B)z_{00}-(1-\frac{1}{2}{T}_B)z_{10}}{(1-\frac{3}{2}{T}_B)}$ (公式2)

边框外的第二层虚拟点的迭代方程：

z_-20＝z₂₀+α²(z₁₁+z_1-1-z_-11-z_-1-1)-2(1+α²)(z₁₀-z_-10)(公式3)

边框外的位于角落的虚拟点的迭代方程：

z_-1-1＝z_1-1+z_-11-z₁₁(公式4)

(4)靠近断层的虚拟点的迭代方程

到目前为止，所有的网格点以及虚拟点都有对应的迭代方程，将断层数据加入进来后，问题变得异常复杂。一条断层线是一个二维的切断文件，它定义了一条在网格化插值时负责阻断信息流的线段，破坏了断层两边的连续性，使得断层的两边不连续。当网格化插值一个数据集时，当计算位于断层线一侧网格点的值时，位于断层线另外一侧的数据不能被直接使用。

如图7所示的黑色矩形所处的网格点在断层的附近，故它在使用公式1时发现在断层左侧的网格点不能使用，因为断层破坏了它们之间的连续性。为此，引入与边框外的虚拟点类似的五个断层附近的虚拟点(这些虚拟点是相对于黑色矩形所处的网格点)，断层附近的网格点的迭代方程可以完全按照公式1进行。

虚拟点迭代计算方法如下：如图8所示，黑色矩形所处的网格点为断层附近的虚拟点，虚拟点所属的每个网格点都有对应的迭代方程，那么我们将这些所有网格点的迭代方程加权平均，就得到了此虚拟点的迭代方程了。这样做目的是保证断层附近的网格点不会出现极值情况，防止随着迭代次数的增加，插值会越来越大(或越来越小)，不会出现收敛的条件。

(5)受原始数据点约束的网格点的迭代方程

如图9所示的受原始数据点约束的网格点示意图，中心网格点(0，0)受原始数据点E约束，它的迭代方程不能简单的使用公式1，原始数据点必须要加入到迭代方程，这样才能更加逼近真实。关于中心网格点使用二次泰勒展开式来近似，同时使用五个不同的点(A，B，C，D四个网格点，加上原始数据点E)。图中的E点位于中心网格点的第一个象限区域中，故选取了对应的五个点；如果E点位于其它的象限中，那么其余四个网格点的选取就类似于旋转得到。设五个不同的实数，b_k(其中k＝1，2，3，4，5)，设 根据公式5所列的方程组，可以求解b_k的值。

$\left[\begin{array}{ccccc}-1& 1& 0& 1& \mathrm{\ξ}\\ 1& 0& -1& -1& \mathrm{\α\η}\\ 1& 1& 0& 1& {\mathrm{\ξ}}^2\\ -1& 0& 0& -1& \mathrm{\ξ\α\η}\\ 1& 0& 1& 1& {\mathrm{\α}}^2{\mathrm{\η}}^2\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\\ b_4\\ b_5\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 2\\ 0\\ 2{\mathrm{\α}}^2\end{array}\right]$ (公式5)

按照差分方程，受原始数据点约束的网格点的迭代方程如下：

$z_{00}=\frac{\left\{(1-T_I)\left\{\begin{array}{c}{z}_{20}+z_{-20}+{\mathrm{\α}}^4(z_{02}+z_{0-2})\\ +2{\mathrm{\α}}^2(z_{11}+z_{1-1}+z_{-11}+z_{-1-1})\\ -2(1+{\mathrm{\α}}^2)[z_{10}+z_{-10}+{\mathrm{\α}}^2(z_{01}+z_{0-1})+\mathrm{\Σ}(b_k\·z_k)]\\ -T_I\mathrm{\Σ}(b_k\·z_k)\end{array}\right\}\right\}}{\{T_I\mathrm{\Σ}{b}_k-2(1-T_I)[(1+{\mathrm{\α}}^4)-(1+{\mathrm{\α}}^2)\mathrm{\Σ}{b}_k\left]\right\}}$ (公式6)

最后，如果E点直接精确地落入到网格点上，则公式5中的矩阵是奇异的，那么公式6是无法直接使用的。但是，这种情况下z₀₀可以简单地用z_E替换。

(6)断层端点附近的网格点的迭代方程

如图10和图11所示的断层端点附近的网格点示意图，灰色圆圈表示的网格点(或是虚拟点)位于断层的端点附近，这些网格点不能简单地认为被断层完全隔离，它们被阻断，但是还有一部分信息交流，也就是说断层端点附近的网格点之间有一定连通性。如何通过迭代方程的建立来描述这种连通性呢，使得问题进入更深层次的讨论。本发明提出的方案是，虚拟点的迭代方程的建立既要和所属的网格点有关，也要和其连通点有关。

如图11所示，A点表示虚拟点，它所属于B，C，D网格点。但是虚拟点A与网格点E有连通关系，故虚拟点A的迭代方程不仅仅由B，C，D的迭代方程密切相关，而且还要加入网格点E的影响。即虚拟点A由B，C，D，E点迭代方程加权平均得到。这样既考虑了断层的非连续性，又考虑了部分连通性。

(7)设置网格点的初始值

在进行迭代求解时，需要设置网格点的初始值。那么设置初始值的原则是尽量的与真实的值接近，同时执行的效率高。一般设置初始值的方法是反距离加权法，但是由于断层线的加入，如果每次反比例加权部检查断层是否阻断，这样在实际情况中的应用极大地影响了效率。

本发明提出“以种子点扩散法”来设置初始值，其基本思路如下：以受原始数据点约束的网格点为种子，向上、下、左、右，以及左上、右上、左下、右下，八个方向扩散，在扩散过程中需要查看对应的边是否为断层边，如果为断层边，则说明被断层阻断，朝此方向的扩散结束。循环每个种子点的扩散，那么每个网格点都会得到一个关于种子点和扩散次数的列表，根据这个列表通过反距离加权法可以计算初始值。这种设置网格点初始值的方法比传统的方法在效率上至少提高10倍。虚拟点也需要初始值，将虚拟点的初始值设置为它所属的所有网格点的加权平均值。

(8)迭代求解

到目前为止，所有的网格点和虚拟点都有对应的迭代方程和初始值。迭代求解采用高斯-赛德尔方法，同时设置最大迭代次数和迭代精度来控制迭代的进行。当达到一定的迭代精度后，再继续迭代意义不大，故选择停止迭代。

第五步，等值线追踪

等值线图包含多条高程值互不相同的等值线，其中每条等值线是由一系列高程值相同的等值点连接而成的。

等值线具有三条基本性质：等值线互不相交；指定高程值的等值线数量可能不止一条；等值线分为闭合等值线和非闭合等值线。

经过网格化插值计算后，每个网格点都有对应的高程值，则整个研究区域形成了一个数据表面，这为等值线的追踪提供了数据支持。在三角形列表中，每个三角形有三条边，并且每条边上都包含有断层信息，边框信息等等。本发明在等值线追踪中采用的基本原理是基于三角形的。

下面以追踪一条高程值为elevation的等值线为例予以说明：

指定边(其两个端点的高程值分别是z₁，z₂)与高程值为elevation的等值线是否相交的判断方法如下：

(elevation-z₁)(elevation-z₂)＜0(公式7)

如果公式7成立，则说明等值线通过该边，接着可以在此边上线性插值获得等值点的具体位置；否则不通过。

遍历每条边，根据公式7判断，如果发现该等值线通过该边时，那么将此边设置为追踪方法中的初始边。下面各个图中，用加粗表示的边称为初始边。接着需要以该初始边为起点，因为一条边有两个三角形与它关联，所以需要从两个方向分别追踪剩余的等值点。等值线首先通过加粗的初始边，形成第一个等值点0点。

对于每个三角形，如果等值线能够从其中一边进入，那么它必然会从其它两边中选择出去。在追踪过程中，需要为每条边设置一个标记，记录此边是否在当前等值线追踪中已经使用过，用于识别封闭等值线的形成。由此寻找下一条边，下一个三角形，依次寻找1点、2点、....、6点。

有三种情况具备追踪结束条件：(1)等值线形成环，形成封闭的等值线。(2)等值线停靠到边框线。(3)等值线停靠到断层线。

下面分别根据以上三种情况依次讨论：

第一种情况，如图12所示的等值线循环追踪示意图，等值线在循环追踪过程中，发现当前正在追踪的边重新回到初始边，则循环立即结束，说明封闭的等值线已经形成。

第二种情况，如图13所示等值线停靠到边框线，等值线在循环追踪过程中，发现当前正在追踪的边是边框边，如图中点6，点6点是在边框边上线性插值形成的，此时循环也应该立即结束，则说明等值线的一个方向的追踪已经结束，那么这条等值线必然是非封闭的等值线。

如图14所示的等值线停靠到断层线的追踪示意图，等值线在追踪过程中，发现当前预通过的边为断层边(如边AB)。此时当前边的两个端点的值不能简单地分别取A，B两个网格点，因为A，B两个网格点由于断层的阻断，它们的网格点的高程值没有任何关系。根据实际情况，等值线来自于0点，1点，故A点应该选择真实的网格点，而B点应该选择A点的虚拟点。在AB线段上根据线性插值，计算出下一个等值点为2点。

同理，等值线预通过A，C两点组成的边，此时A点取真实的网格点，而C点取A点的虚拟点，经过线性插值得到4点，但是发现4点和2点连成的线段与断层线相交于3点，那么最后一个等值点已经找到，它就是3点。追踪到此结束。

同时注意：正在探测的等值点与上一个等值点有可能与多条断层线相交，这时候需要取离上一个等值点最近的那个相交点作为等值线的结束点。

此等值线正好停靠在断层线上，由于采用了虚拟点技术，这种停靠方法更加自然，避免了其它方法中停靠不精确的情况。

第六步：等值线填充区域的搜索

等值线图在平面内无法反映高程值整体变化趋势，故引入等值线区域的填充，它能够直观地反映趋势变化。等值线在边框线、断层线的约束下会形成各种闭合的区域，这些区域可以填充不同的渐变颜色或纹理给用户直观的感受。

每一个填充区域部是封闭的环结构，故在填充区域的搜索中，我们也建立对应的环结构。

如图15所示的基于边框搜索线和断层搜索线形成的环结构示意图，所有的内部边框线按照逆时针方向形成一个环，称为边框搜索线环。同时其它的断层搜索线部是相对于边框搜索线环而言是相反的。

建立两类不同的区域搜索线：边框搜索线，断层搜索线。如图15所示：直线AB，BC，CD，DA都是边框搜索线；曲线EFGHI，JFKLM为断层搜索线。同时建立搜索点列表，搜索点的类型包含五类：边框线的端点(如A，B，C，D)，边框线与等值线的交点，断层线的折点(如G，H，K，L)，断层线的端点(如E，I，J，M)，断层线与等值线的交点，断层线与断层线的交点(如F)。

如图16所示的由封闭的等值线形成的填充区域示意图，由封闭的等值线可以直接形成区域1，区域2，它们也是最简单的填充区域。但是我们发现区域1包含在区域2中，在成图的过程中需要利用这样的包含关系。在绘图时，需要先填充区域2的颜色或纹理，再填充区域1的颜色或纹理。

如图17所示的由非封闭的等形成的填充区域示意图，由非封闭的等值线围成的填充区域如区域3，区域4，区域5，区域6。这些区域的搜索相对比较复杂。需要遍历所有的边框搜索线，断层搜索线。现在以区域3的搜索过程为例予以说明。发现边框搜索线上的边AB没有搜索，则此AB边必然对应了一个填充区域，由此以它为起始边，开始区域搜索。检查B点是边框线和等值线的交点，则顺着等值线到达C点。又C点是等值线和断层线的交点，则顺着断层搜索线到达D，E，F点。类似地，发现F点来自断层与另外一条等值线，则顺着等值线来到A点。同时，发现A点是这次搜索的起始点，故区域搜索结束。那么形成区域3所有控制点就找到了，完成了区域3的搜索。

第七步：等值线数据的格式化

网格化插值、等值线的追踪和区域的搜索都需要复杂的计算，需要耗费大量的内存和时间，因此等值线生成后需要把数据格式化后保存下来，以便入库或保存到本地文件中。这样就避免了每次打开等值线图时都重新计算。

通过本发明的基于最小曲率法插值的含断层的等值线图自动生成方法，能够实现如图18所示的等值线图最终生成的总体效果图。

本发明技术方案改进的效果如下：

1.等值线图能够快速生成；

2.等值线追踪算法更加简单、容易实现；

3.等值线的追踪算法基于三角形，在三角形中追踪等值点可以避免在四边形追踪中会出现选择困难的问题；

4.等值线能够精确地停靠到断层线，如图19所示的等值线的断层停靠效果图，黑色椭圆所指的区域说明等值线在追踪过程中停靠到断层线上时自然且效果美观，由于采用了虚拟点技术，避免了因为停靠不精确而出现等值线在断层线附近的弯折现象；

5.等值区域的纹理填充，如图20所示的等值区域的纹理填充效果图，用平铺的纹理填充等值区域，其它的区域依然用渐变颜色填充。因为等值区域的搜索基于边框或断层的搜索线，能够找到完整的等值区域，避免了其它方法中无法进行纹理填充的问题。

以上通过具体的和优选的实施例详细的描述了本发明，但本领域技术人员应该明白，本发明并不局限于以上所述实施例，凡在本发明的精神和原则之内，所作的任何修改、等同替换等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于最小曲率法插值的含断层的等值线图自动生成方法，其特征在于，包括：

步骤1，网格化，确定研究区域，并在所述研究区域内建立网格；

步骤2，三角化，将每个所述网格以其中一条对角线划分为两个三角形，并将整个所述研究区域建立一个互不交叉的三角形列表；

步骤3，数据的预处理，建立点列表、边列表和三角形列表，同时记录边类型；

步骤4，网格点的插值，采用最小曲率法插值方法对每个所述网格点插值；

步骤5，等值线追踪，采用虚拟点技术进行所述等值线追踪；

步骤6，填充区域的搜索，建立基于边框搜索线和断层搜索线形成的环结构，按照逆时针方向遍历搜索线寻找所述填充区域；

步骤7，等值线数据的格式化，所述等值线数据生成后，将所述等值线数据和所述填充区域的数据格式化后保存；

步骤3中所述的边类型包括边框边、内部边、断层边和虚拟边，所述边框边是指边位于边框上；所述断层边是指边有断层通过；所述虚拟边是指边本身不存在，随着等值线的追踪动态创建；所述内部边是指除所述边框边、所述断层边和所述虚拟边的边；

步骤4中所述网格点插值主要包括：

(1)所述网格点的迭代方程的建立；

(2)所述网格点的初始值设置；

(3)所述网格点的迭代方程的迭代求解；

不同所述网格点的迭代方程包括：分别是远离边框的网格点迭代方程，如公式a；边框外的虚拟点的迭代方程；靠近断层的虚拟点的迭代方程；受原始数据点约束的网格点的迭代方程，如公式d；和断层端点附近的网格点的迭代方程；

所述公式a，如下：

公式(a)

其中，α为网格的横纵比Δx为网格矩形的长，Δy为网格矩形的宽，T_I为内部张力，默认为0；

所述边框外的虚拟点的迭代方程，包括：框外的第一层虚拟点的迭代方程，如公式b-1；边框外的第二层虚拟点的迭代方程，如公式b-2；边框外的位于角落的虚拟点的迭代方程，如公式b-3：

公式(b-1)

z_-20＝z₂₀+α²(z₁₁+z_1-1-z_-11-z_-1-1)-2(1+α²)(z₁₀-z_-10)公式(b-2)

z_-1-1＝z_1-1+z_-11-z₁₁公式(b-3)

所述靠近断层的虚拟点的迭代方程:所述虚拟点所属的每个网格点都有对应的迭代方程加权平均，得到所述虚拟点的迭代方程；

公式(d)

所述断层端点附近的网格点的迭代方程：虚拟点A的迭代方程由所述虚拟点A所属的网格点的迭代方程和与所述虚拟点A有连通关系的网格点的迭代方程加权平均得到；

其中，z₀₁，z_0-1，z₁₀，z_-10表示目标网格点的四个内环点的高程值，z₁₁，z_-11，z_1-1，z_-1-1表示目标网格点的四个中环点的高程值，z₀₂，z_0-2，z₂₀，z_-20表示目标网格点的四个外环点的高程值；

T_B表示边框张力，取值范围是0到1；

z₁、z₂、z₃和z₄表示相对于原始数据点和目标网格点的临近网格点的高程值，z₅表示原始数据点的高程值；

b₁、b₂、b₃、b₄表示临近网格点相对于目标网格点的膨胀因子，b₅表示原始数据点相对于目标网格点的膨胀因子；

在建立所述断层端点附近的网格点的迭代方程时，通过周围所述网格点加权平均来设置；

建立边框附近的所述网格点的迭代方程和断层附近的所述网格点迭代方程时，当出现邻近网格点不足的情况时，附加两排虚拟网格点，通过所述虚拟网格点所属的所述网格点的迭代加权平均来建立所述虚拟网格点的迭代方程；

所述网格点的初始值的设置方法如下：以受所述原始数据点约束的所述网格点为种子，向上、下、左、右以及左上、右上、左下、右下八个方向扩散，在扩散过程中查看对应的边是否为断层边，如果为断层边，则说明被断层阻断，朝此方向的扩散结束；循环每个种子点的扩散，每个所述网格点都会得到一个关于所述种子点和扩散次数的列表，根据所述列表，通过反距离加权法计算所述网格点的初始值；

所述等值线追踪采用基于三角形的基本原理，同时采用所述虚拟点技术，拓展断层单侧数据的连续性，在所述等值线追踪过程中，时刻检查所述等值线是否和断层相交，并选择所述等值线和断层相交的点做为停靠点。

2.根据权利要求1所述的基于最小曲率法插值的含断层的等值线图自动生成方法，其特征在于，所述研究区域为包含原始数据点和断层线的最大矩形，所述原始数据点包含的属性有横坐标x、纵坐标y、高程值z。

3.根据权利要求1所述的基于最小曲率法插值的含断层的等值线图自动生成方法，其特征在于，建立所述网格点的迭代方程时，首先建立以所述网格点为中心的中心区域，在所述中心区域内，通过将多个所述原始数据点加权平均，使得每个所述网格点的矩形区域至多有一个原始数据点与之对应。