CN104182571A - 基于Delaunay和GPU的Kriging插值方法 - Google Patents

Abstract

本发明公开了一种基于Delaunay和GPU的Kriging插值方法。其包括以下步骤：构建三维Delaunay三角网、建立空间索引、定位插值点所在四面体和实现Kriging插值运算。本发明的有益效果是：解决了Kriging插值过程中邻域搜索和插值运算速度慢的问题，具有效率高的优点；同时，当采样点较多时，本发明具有更好的适应性。

Description

基于Delaunay和GPU的Kriging插值方法

技术领域

本发明属于插值方法技术领域，尤其涉及一种基于Delaunay和GPU的Kriging插值方法。

背景技术

数据可视化技术指的是运用计算机图形学和图像处理技术等技术，将数据转换为图形或图像在屏幕上显示出来，并进行交互式处理的相关理论、方法和技术。它涉及到的领域广泛，如计算机辅助设计、计算机视觉、人机交互技术、计算机图形学、图像处理等多个领域。基于三维空间的数据可视化方法在计算机图形技术发展的影响下也越来越普及。Kriging插值算法是实现数据可视化的行之有效的方法和工具，在地理信息、气象预报、油气勘探、图形图像处理及工程可视化等多个领域发挥着举足轻重的作用。Kriging插值算法是一种线性无偏最优插值方法，主要用来处理空间现象，它由矿业工程师Krige首次提出。克里金插值的表达形式很简单，即对于未知点处的估值，采用的是对其周围可用点进行线性组合的形式。克里金方法的主要任务是：定量的估计每个已知点的属性值对未知点属性值的影响。该影响的大小用权重系数来表示，最后利用这些权重系数进行加权平均得到线性、无偏、和最小方差的未知点的属性值。近几十年来，Kriging方法在采矿学、环境科学、医学、土壤学等研究领域都有着极其广泛的应用。点集的三角剖分(Triangulation)，对数值分析(比如有限元分析)以及图形学来说，都是极为重要的一项预处理技术。尤其是Delaunay三角剖分，由于其独特性，关于点集的很多种几何图都和Delaunay三角剖分相关，如Voronoi图，EMST树，Gabriel图等。Delaunay三角剖分有最大化最小角，“最接近于规则化的“的三角网和唯一性(任意四点不能共圆)两个特点。逐点插入的Lawson算法是Lawson在1977年提出的，该算法思路简单，易于编程实现。基本原理为：首先建立一个大的三角形或多边形，把所有数据点包围起来，向其中插入一点，该点与包含它的三角形三个顶点相连，形成三个新的三角形，然后逐个对它们进行空外接圆检测，同时用Lawson设计的局部优化过程LOP进行优化，即通过交换对角线的方法来保证所形成的三角网为Delaunay三角网。CUDA是一种由NVIDIA推出的通用并行计算架构，该架构使GPU能够解决复杂的计算问题。它包含了CUDA指令集架构(ISA)以及GPU内部的并行计算引擎。支持CUDA的GPU销量已逾1亿，数以千计的软件开发人员正在使用免费的CUDA软件开发工具来解决各种专业以及家用应用程序中的问题。这些应用程序从视频与音频处理和物理效果模拟到石油天然气勘探、产品设计、医学成像以及科学研究，涵盖了各个领域。GPU已经不再局限于3D图形处理了，GPU通用计算技术发展已经引起业界不少的关注，事实也证明在浮点运算、并行计算等部分计算方面，GPU可以提供数十倍乃至于上百倍于CPU的性能。在Windows7中，CPU与GPU组成了协同处理环境。CPU运算非常复杂的序列代码，而GPU则运行大规模并行应用程序。然而，随着工程应用的不断深入，为了达到更精确的数据可视化效果，Kriging插值所需的插值点数量是庞大的，这必然导致插值效率下降。目前，由于GPU技术的快速发展，已有多种利用GPU来提高Kriging插值速度的算法被提出，但其并不适用于所有情况，尤其当采样点较多时。

发明内容

为了解决以上问题，本发明提出了一种基于Delaunay和GPU的Kriging插值方法。

本发明的技术方案是：一种基于Delaunay和GPU的Kriging插值方法，包括以下步骤：

S1.采用Delaunay三角剖分方法，将采样点构建为三维Delaunay三角网；

S2.利用S1中构建的三维Delaunay三角网，建立索引，定位包含待插值点的四面体；

S3.通过CUDA在GPU上进行并行运算，实现Kriging插值运算，具体包括以下步骤：

S3.1.将Kriging方程组转换为矩阵形式表示为：

[K]·[λ]＝[M]

其中，K矩阵由每两个已知点之间的变差函数值组成，M矩阵由当前待插值点与所有与知点之间的变差函数值组成，λ矩阵为权重系数矩阵；

S3.2.将数据从主机端拷贝到设备端，包括已知点坐标、属性值、领域点坐标及待插值点坐标信息；

S3.3.通过kernel函数为每一个待插值点分配一个线程，通过并行矩阵运算得到待插值点属性值；

S3.4.将S3.3中得到的属性值从设备端拷贝回主机端，实现Kriging插值运算。

进一步地，所述步骤S2中建立索引包括建立八叉树索引和建立空间索引。

进一步地，所述建立八叉树索引的算法具体包括以下步骤：

S2.1.设定阀值N和插值区域S，并将区域内的所有四面体放入根节点中；

S2.2.若根节点内的四面体个数小于N，则完成索引建立；

S2.3.若根节点内的四面体个数大于N，则将根节点八等份为八个子区域，生成八个子节点；

S2.4.依次将父节点中的四面体移入它包含的子节点中，并记录各个子节点

区域内的四面体个数；

S2.5.依次检查八个子节点中四面体个数，若子节点内的四面体个数大于N，则重复步骤S2.3；

S2.6.若子节点内的四面体个数小于N，则完成索引建立。

进一步地，所述建立空间索引的算法具体包括以下步骤：

S2.7.设定插值区域S，从区域内选取一个四面体，计算该四面体四个顶点的最大Y值、最小Y值、最大Z值和最小Z值，并计算得到Z方向的索引号范围；

S2.8.将S2.7中得到的Z值按照从小到大的顺序固定，并依次对该四面体进行切割，计算该Z平面与四面体的交点个数；

S2.9.若S2.8中计算得到的交点个数为1，则完成索引建立；

S2.10.若S2.8中计算得到的交点个数为3或4，则将交点按照Y值大小排序，并将交点构成三角形或四边形；

S2.11.根据S2.7中最大Y值和最小Y值计算得到Y方向的索引号范围，将Y值按照从小到大的顺序固定；

S2.12.依次对S2.10中的三角形或四边形进行切割，计算该平行于X轴的直线与三角形或四边形的交点(X₁，X₂)，则落在X₁与X₂范围内的点即属于该四面体，完成索引建立。

本发明的有益效果是：本发明的基于Delaunay和GPU的Kriging插值方法通过为三角网建立空间索引的方法和基于CUDA的并行插值策略，解决了Kriging插值过程中邻域搜索和插值运算速度慢的问题，具有效率高的优点；同时，当采样点较多时，本发明具有更好的适应性。

附图说明

图1是本发明的基于Delaunay和GPU的Kriging插值方法流程示意图。

图2是本发明的三维Delaunay三角网结构示意图。

图3是本发明的建立八叉树索引算法流程示意图。

图4是本发明建立的八叉树索引结构示意图。

图5是本发明的建立空间索引算法流程示意图。

图6是本发明的并行矩阵运算的流程示意图。

图7是本发明的并行矩阵运算实现代码图。

图8是本发明的三维Kriging插值并行实现效果图。

图9是现有方法与本发明方法的误差对比图。

图10是本发明实施例中的对比图。

具体实施方式

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

如图1所示，为本发明的基于Delaunay和GPU的Kriging插值方法流程示意图。本发明结合了计算几何领域Delaunay三角网技术和Kriging插值算法，通过对采样点构建三维Delaunay三角网，进而对三维三角网建立空间索引，以实现快速搜索邻域点，再利用GPU强大的并行计算能力，借助CUDA平台实现高效率的插值运算，克服了现有技术中将所有采样点作为已知点来进行插值计算时效率不高或无法实现的问题。本发明的基于Delaunay和GPU的Kriging插值方法包括以下步骤：

S1.采用Delaunay三角剖分方法，将采样点构建为三维Delaunay三角网。

如图2所示，为50个点构成的Delaunay三角网。Delaunay三角剖分方法主要用于将空间中散乱的采样点规则化；通过对其构建三维Delaunay三角网，可以合理有效的对这些采样点的空间位置进行记录和组织。

S2.利用S1中构建的三维Delaunay三角网，建立索引，定位包含待插值点的四面体。

在Kriging插值过程中，如果对每一个待插值点都要遍历所有的四面体，将会严重影响插值效率。因此本发明通过步骤S1构建三维Delaunay四面体网格，并对采样点的空间位置进行记录和组织。这样我们可以通过建立相应的索引机制来实现点的快速定位。对于索引机制的选取，需要依据插值点个数和四体面个数来决定。当待插值点个数较少而四面体个数较多时，我们选择八叉树索引算法。因为基于八叉树索引的邻域搜索算法需要遍历所有待插值点，当待插值点较少时，运算速度较快。如图3所示，为本发明的建立八叉树索引算法流程示意图。建立八叉树索引的算法具体包括以下步骤：

S2.1.设定阀值N和插值区域S，并将区域内的所有四面体放入根节点中；

S2.2.若根节点内的四面体个数小于N，则完成索引建立；

S2.3.若根节点内的四面体个数大于N，则将根节点八等份为八个子区域，生成八个子节点；

S2.4.依次将父节点中的四面体移入它包含的子节点中，并记录各个子节点区域内的四面体个数；

S2.5.依次检查八个子节点中四面体个数，若子节点内的四面体个数大于N，则重复步骤S2.3；

S2.6.若子节点内的四面体个数小于N，则完成索引建立。

如图4所示，为本发明建立的八叉树索引结构示意图。在构建八叉树的过程中，需要注意：八叉树的建立，实际上是将数据所属区域按四面体个数不超过N的原则，不断对其八等分的过程。因此，八叉树的深度与阈值N有很大的关系。N越大，八叉树的深度越小，构建索引时间越短，但节点内搜索四面体效率越低；反之N越小，八叉树深度越大，索引构建时间越长，但节点内搜索四面体效率提高。为了合理选取N的值来平衡八叉树的构建时间和检索时间，我们根据四面体个数和待插值点个数的多少来选取N的值，待插值点数多时N选小一些，以减少每个采样点检索四面体的个数，四面体个数多时N应选大一些，避免八叉树深度过大，构建索引时间太长。

当待插值点个数较多而四面体个数较少时，由于基于八叉树索引的邻域搜索算法需要遍历所有待插值点，运算速度会受到影响。因此，我们选择基于待插值点的空间索引算法。当待插值点较多时，由于每一个待插值点只属于一个四面体，我们依次取出一个四面体，找出所有落在该四面体内的待插值点，遍历完所有的四面体，就完成了所有待插值点的定位。如图5所示，为本发明的基于待插值点建立空间索引算法流程示意图。基于待插值点建立空间索引的算法具体包括以下步骤：

S2.7.设定插值区域S，从区域内选取一个四面体，计算该四面体四个顶点的最大Y值、最小Y值、最大Z值和最小Z值，并计算得到Z方向的索引号范围；

S2.8.将S2.7中得到的Z值按照从小到大的顺序固定，并依次对该四面体进行切割，计算该Z平面与四面体的交点个数；

S2.9.若S2.8中计算得到的交点个数为1，则完成索引建立；

S2.10.若S2.8中计算得到的交点个数为3或4，则将交点按照Y值大小排序，并将交点构成三角形或四边形；

S2.11.根据S2.7中最大Y值和最小Y值计算得到Y方向的索引号范围，将Y值按照从小到大的顺序固定；

S2.12.依次对S2.10中的三角形或四边形进行切割，计算该平行于X轴的直线与三角形或四边形的交点(X₁，X₂)，则落在X₁与X₂范围内的点即属于该四面体，完成索引建立。

我们通过对三维Delaunay三角网建立索引机制，可以排除大量与待插值点无关的四面体，将要定位的目标锁定在少量的四面体内，缩小搜索的范围。

对于每一个待插值点相同的计算过程，我们在应用程序中将高密度数据进行并行计算的部分做成一个称作kernel的函数在GPU设备上执行，通过设计kernel中的线程来完成通用计算的GPU实现。如图6所示，为本发明的并行矩阵运算的流程示意图。

S3.通过CUDA在GPU上进行并行运算，实现Kriging插值运算，具体包括以下步骤：

S3.1.在构建好空间索引后，将Kriging方程组转换为矩阵的形式，由λ＝[K]⁺[M]计算得到：

[K]·[λ]＝[M]

其中，K矩阵由每两个已知点之间的变差函数值组成，M矩阵由当前待插值点与所有与知点之间的变差函数值组成，λ矩阵为权重系数矩阵；

S3.2.将数据从主机端拷贝到设备端，包括已知点坐标、属性值、领域点坐标及待插值点坐标信息；

S3.3.通过kernel函数为每一个待插值点分配一个线程，通过并行矩阵运算得到待插值点属性值；

如图5所示，为本发明的并行矩阵运算实现代码图。N表示已知点个数，M表示待插值点个数。在kernel中，为每一个待插值点分配一个线程，在线程内通过一系列的矩阵运算得到待插值点属性值。

S3.4.将S3.3中得到的属性值从设备端拷贝回主机端，实现Kriging插值运算。

如图7所示，为本发明的并行矩阵运算实现代码图。这里，我们给出并行实现策略。首先，将数据从主机端拷贝到设备端；然后主机启动kernel，在kernel中，为每一个待插值点分配一个线程，在线程内通过一系列的矩阵运算得到待插值点属性值；最后，将数据从设备端拷贝回主机端。

基于以上提出的并行插值算法，本专利实现了对采样得到的实测数据的三维Kriging插值运算。该实测数据包括坐标值、深度值、速度值等。速度可以看成是一种区域化变量，坐标和深度值组成已知点的空间坐标值，已知点的属性值即为速度值。如下表1所示为已知点数据信息。

坐标	最小值	最大值
			x方向	18486665.5	18659409
y方向	3305982	3400572
			z方向	909.06	1875

如图8所示，为本发明的三维Kriging插值并行实现效果图。已知点个数为10000，插值点个数为8120601。如图9所示，为利用所有已知点插值与本发明插值结果的误差对比图。经计算得到最大相对误差为0.83％。

下面我们对本发明提出的基于Delaunay和GPU的Kriging插值方法的插值效率进行了实验。实验实现的环境为：Windows7，Intel Core i33.4GHz CPU，4.0G RAM，NVIDIA GeforceGT 630。我们分别用现有的GPU插值算法与本发明的方法对以上提到的实测数据进行Kriging插值运算。在现有的GPU插值方法中，矩阵的相乘运算部分是用CUBLAS库中提供的cublasDgemm函数来实现的，而矩阵赋值、求逆等其它部分则在CPU上实现。如图10所示，为当已知点固定在103，待插值点个数不同时，两种插值算法的时间对比及时间差。通过图10可知，本发明提出的基于Delaunay和GPU的Kriging插值方法大大减少了插值时间，尤其随着插值点数的增加，效率优势更加明显。如下表2所示，为当待插值点固定在33201，随着已知点个数的增加，两种插值方法的时间对比。

其中，GPU1 time代表基于已有的GPU插值算法的运算时间，GPU2 time代表基于本文提出的插值算法的运算时间，N/A表示运算时间过长或者已知点个数太大而无法进行运算。由表2可知，本发明提出的基于Delaunay和GPU的Kriging插值方法明显优于现有的基于GPU的插值方法；并且，当已知点在到时一定的值时，现有的基于GPU的插值方法是不可行的，因此本发明具有更好的适应性。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (4)

1.一种基于Delaunay和GPU的Kriging插值方法，其特征在于，包括以下步骤：

S1.采用Delaunay三角剖分方法，将采样点构建为三维Delaunay三角网；

S2.利用S1中构建的三维Delaunay三角网，建立索引，定位包含待插值点的四面体；

S3.通过CUDA在GPU上进行并行运算，实现Kriging插值运算，具体包括以下步骤：

S3.1.将Kriging方程组转换为矩阵形式表示为：

[K]·[λ]＝[M]

其中，K矩阵由每两个已知点之间的变差函数值组成，M矩阵由当前待插值点与所有与知点之间的变差函数值组成，λ矩阵为权重系数矩阵；

S3.2.将数据从主机端拷贝到设备端，包括已知点坐标、属性值、领域点坐标及待插值点坐标信息；

S3.3.通过kernel函数为每一个待插值点分配一个线程，通过并行矩阵运算得到待插值点属性值；

S3.4.将S3.3中得到的属性值从设备端拷贝回主机端，实现Kriging插值运算。

2.如权利要求1所述的基于Delaunay和GPU的Kriging插值方法，其特征在于：所述步骤S2中建立索引包括建立八叉树索引和建立空间索引。

3.如权利要求2所述的基于Delaunay和GPU的Kriging插值方法，其特征在于：所述建立八叉树索引的算法具体包括以下步骤：

S2.1.设定阀值N和插值区域S，并将区域内的所有四面体放入根节点中；

S2.2.若根节点内的四面体个数小于N，则完成索引建立；

S2.3.若根节点内的四面体个数大于N，则将根节点八等份为八个子区域，生成八个子节点；

S2.4.依次将父节点中的四面体移入它包含的子节点中，并记录各个子节点

区域内的四面体个数；

S2.5.依次检查八个子节点中四面体个数，若子节点内的四面体个数大于N，则重复步骤S2.3；

S2.6.若子节点内的四面体个数小于N，则完成索引建立。

4.如权利要求2所述的基于Delaunay和GPU的Kriging插值方法，其特征在于：所述建立空间索引的算法具体包括以下步骤：

S2.7.设定插值区域S，从区域内选取一个四面体，计算该四面体四个顶点的最大Y值、最小Y值、最大Z值和最小Z值，并计算得到Z方向的索引号范围；

S2.8.将S2.7中得到的Z值按照从小到大的顺序固定，并依次对该四面体进行切割，计算该Z平面与四面体的交点个数；

S2.9.若S2.8中计算得到的交点个数为1，则完成索引建立；

S2.10.若S2.8中计算得到的交点个数为3或4，则将交点按照Y值大小排序，并将交点构成三角形或四边形；

S2.11.根据S2.7中最大Y值和最小Y值计算得到Y方向的索引号范围，将Y值按照从小到大的顺序固定；

S2.12.依次对S2.10中的三角形或四边形进行切割，计算该平行于X轴的直线与三角形或四边形的交点(X₁，X₂)，则落在X₁与X₂范围内的点即属于该四面体，完成索引建立。