CN105653501B - 一种加速克里金插值的方法 - Google Patents

Abstract

一种加速克里金插值的方法，涉及地学领域。本发明为解决现有加速克里金插值的方法仅从计算平台的计算资源的角度来加速，克里金插值中仍存在大量的冗余计算和存储的问题。一种加速克里金插值的方法，主要按以下步骤进行：搜索待估未知点x周围n个已知观测点；并标识具有相同邻近观测点的待估未知点，得到父节点；标记为父节点的GPU线程，及其拟合模型构造变异函数矩阵K′；标记为父节点的GPU线程，求K′的逆矩阵K′^‑1；GPU内所有线程构造自己的矩阵M′；GPU内所有线程解方程序组求出λ′；GPU内所有线程根据系数λ_i估计待测值x。本发明适用于地学领域。

Description

一种加速克里金插值的方法

技术领域

本发明涉及地学领域，尤其涉及一种加速克里金插值的方法。

背景技术

克里金插值算法也称为空间自协方差最佳插值法，它最早由南非矿业工程师D.G.Krige提出，因此也以他的名字命名该算法。克里金插值算法广泛用于地学、环境科学等领域。它是根据未知点周围一定范围内观测点的值来估计未知点取值的一种方法。该方法是一种线性无偏方法，内插结果准确率高。然而在实际应用中，待插值地区面往往过大，插值密度也要求很高，这就造成了克里金插值速度很慢，严重限制了它的应用。随着并行计算技术的发展，基于这些新技术加速克里金算法成为一个研究热点。目前，克里金加速算法主要分为以下几个方面：

一、基于集群技术

集群由大量计算机通过高速网络连接组成，普通的商用集群一般包含上百个节点。集群的特点是相对普通PC(个人电脑)，运算能力比较强，但也存在着能耗高，管理成本高，应用难度高的缺点。基于集群硬件平台来加速克里金插值方法通常使用消息传递接口(MPI)，具体实现方法包括两种：一是将算法中每一个可以并行的部分在集群中以并行方式完成；二是以每个未知点的估计做为一个并行单元在集群上并行执行。前者编程相对容易，因为克里金插值中的可并行部分是矩阵计算，现有大量专业人员编写的可在集群上并行运行的代码库可以使用，但是这种并行方式对集群的利用率不高，同时会产生大量的节点间通讯，对速度有很大的影响；后者能充分利用集群的计算资源，但是需要编写克里金算法中的包括矩阵计算的各种数值算法，因此实现成本高，对编写代码人员的编程能力要求较高。

二、基于多核CPU

由于工艺原因，CPU的核心频率遇到了瓶颈，CPU的主流发展朝多核方向发展，充分利用多核CPU的资源也成为加速克里金的一种主要方法。多核CPU成本较集群低，实现相对容易，节点间通讯速度也较快，但受CPU核心数量的限制，其能够并行执行的线程数有限，因此加速比潜力不如集群或GPU。基于多核CPU加速克里金插值一般使用OpenMP框架，具体的加速方法也分为两类：一是将可并行执行的矩阵计算在多核CPU上并行执行；二是以每个待估的未知点做为基本的并行单元完成克里金插值的并行执行。

三、基于GPU

GPU平台是一种有前景的并行加速平台。GPU相对集群成本低很多，能耗也小很多，而且GPU包含的可并行执行的计算核心数量远高于多核CPU，主流的GPU都拥有几千个计算核心。基于GPU硬件平台的加速技术主要基于两种软件平台：统一计算设备架构(CUDA)和开放运算语言(OpenCL)。具有加速方法同集群和多核CPU一样，分为两类：一是将可并行执行的矩阵计算在多核CPU上并行执行；二是以每个待估的未知点做为基本的并行单元完成克里金插值的并行执行。

四、组合方法

由于一种计算设备可能同时包含上述几种计算平台，因此上述几种加速算法也可以混合应用，以最大程度挖掘计算平台的潜力。

发明内容

克里金插值是地学领域经常用到的一种插值方法，同其它插值方法相比，例如：反距离加权插值法、最近邻点插值法、线性插值三角网法、局部多项式法等，其插值精度较高。但是克里金插值是一种严重消耗计算和内存资源的算法，尤其是在大面积区域中做高密度克金金插值是相当慢的，严重地影响了克里金插值的实际应用价值。本发明针对克金格插值算法提出一种新的解决方案，能极大的提高克里金插值的速度。

本发明为解决现有加速克里金插值的方法仅从计算平台的计算资源的角度来加速，克里金插值中仍存在大量的冗余计算和存储的问题，而提出一种加速克里金插值的方法。

一种加速克里金插值的方法，按以下步骤进行：

步骤一、在CPU端，按已知观测点的横坐标值对观测点从小到大排序；

步骤二、将排好序的数据传送到GPU，并启动GPU端克里金插值程序；

步骤三、在GPU中每个线程执行一个待估点的计算，所有线程首先搜索待估未知点x周围n个已知观测点x_j,j＝1,2,^…,n，作为该未知点的邻近点；

步骤四、同步GPU中同一个workgroup内所有线程，待workgroup内所有的待估未知点的邻近点搜索完成后，执行下一步；

步骤五、所有workgroup内的第一个线程比较所有相邻待估未知点的邻近点，并标识具有相同邻近观测点的待估未知点；

其中，对于具有相同邻近点的待估点中最先出现的点，称为父节点；

步骤六、标记为父节点的GPU线程，根据及其拟合模型构造变异函数矩阵：

在这个步骤中，只有标记为父节点的线程参与了计算并保存结果矩阵，其它线程既没有消耗计算，也没有浪费存储空间；

步骤七、标记为父节点的GPU线程求K′的逆矩阵K′^-1；

步骤八、GPU内所有线程构造自己的矩阵：

步骤九、GPU内所有线程解方程组λ′＝M′·K′^-1，求出λ′，其中，μ是拉格朗日乘子；

步骤十、GPU内所有线程根据系数λ_i估计待测值x

本发明包括以下有益效果：

1、传统方法仅从计算平台的计算资源的角度来加速，本发明方法从根本上解决了克里金插值法的运算速度慢的问题，在克里金插值过程中，选择离待估计最近的一定数量观测点来计算权值系数，相邻的待估点通常有相同的邻近已知观测点，由于待估未知点的取值完全取决于其邻近已知观测点的值，因此，对于具有相同邻近观测点的待估点，只需计算其中一个待估点，由此减少大量的冗余计算和中间计算过程中的矩阵存储量。

2、本发明方法提出一种搜索邻近点的方法，减少了搜索范围，从而提高了搜索速度。

附图说明

图1为具有相同邻近观测点的待估未知点的标识过程示意图。

具体实施方式

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

具体实施方式一、本实施方式所述的一种加速克里金插值的方法，按以下步骤进行：

步骤一、在CPU端，按已知观测点的横坐标值对观测点从小到大排序；

步骤二、将排好序的数据传送到GPU，并启动GPU端克里金插值程序；

步骤三、在GPU中每个线程执行一个待估点的计算，所有线程首先搜索待估未知点x周围n个已知观测点x_j,j＝1,2,^…,n，作为该未知点的邻近点；

步骤四、同步GPU中同一个workgroup内所有线程，待workgroup内所有的待估未知点的邻近点搜索完成后，执行下一步；

步骤五、所有workgroup内的第一个线程比较所有相邻待估未知点的邻近点，并标识具有相同邻近观测点的待估未知点；

其中，对于具有相同邻近点的待估点中最先出现的点，称为父节点；

步骤六、标记为父节点的GPU线程，根据及其拟合模型构造变异函数矩阵：

在这个步骤中，只有标记为父节点的线程参与了计算并保存结果矩阵，其它线程既没有消耗计算，也没有浪费存储空间；

步骤七、标记为父节点的GPU线程求K′的逆矩阵K′^-1；

步骤八、GPU内所有线程构造自己的矩阵：

步骤九、GPU内所有线程解方程序组λ′＝M′·K′^-1，求出λ′，其中，μ是拉格朗日乘子；

步骤十、GPU内所有线程根据系数λ_i估计待测值x

本实施方式包括以下有益效果：

1、传统方法仅从计算平台的计算资源的角度来加速，本实施方式从根本上解决了克里金插值法的运算速度慢的问题，在克里金插值过程中，选择离待估计最近的一定数量观测点来计算权值系数，相邻的待估点通常有相同的邻近已知观测点，由于待估未知点的取值完全取决于其邻近已知观测点的值，因此，对于具有相同邻近观测点的待估点，只需计算其中一个待估点，由此减少大量的冗余计算和中间计算过程中的矩阵存储量。

2、本实施方式提出一种搜索邻近点的方法，减少了搜索范围，从而提高了搜索速度。

具体实施方式二、本实施方式是对具体实施方式一所述的一种加速克里金插值的方法的进一步说明，步骤三中所述搜索过程如下：

1、将待估计点x横坐标值同已经排好序的观测点的横坐标比较，找到其在观测点中的排序，即x_i<x<x_i+1；

2、定义一大小为n的数组用于存储x的n个邻近点，以x点所在位置开始交替向左右两侧搜索，分两种情况处理：1、数组未被填满：将搜索到的点填入数组；2、数组已经被填满：如果搜索到的点满足下述条件1则替换数组内距x最大距离的邻近点，否则继续搜索；如果满足条件2或条件3则完成搜索；

其中，条件1为：x到当前搜索到的观测点的距离小于数组内所有存储的观测点到x的距离；

条件2为：x到当前搜索点的横坐标间的距离大于数组内所有存储的观测点到x的距离；

条件3为：所有观测点搜索完成。

具体实施方式三、本实施方式是对具体实施方式一或二所述的一种加速克里金插值的方法的进一步说明，步骤五中具有相同邻近观测点的待估未知点的标识方法：如图1所示，图中的每一个方格表示一个待估点，A、B…分别表示待估点的邻近点集合，如果相邻的待估点具有同样的邻近点，那么在图中的方格中的字母也相同；

a、分别将每一个待估点同它的左侧的待估点比较(最左侧的待估点除外)，如相同，在图中表示为←；

b、分别将每一个待估点同它的上边的待估点相比较(最上层待估点除外)，如相同，在图中表示为↑；

c、对于具有相同邻近点的待估点中最先出现的点(从左向右，从上向下数)，即图1中b图中标为字母的待估点，称为父节点，并按从左到右，从上到下的顺序标记顺序号，如图1中c图所示；

d、将其它非首次出现的待估点称为首次出现的待估点的子节点，将其填上其父节点的顺序号。

具体实施方式四、本实施方式是对具体实施方式一至三之一所述的一种加速克里金插值的方法的进一步说明，本发明方法适合于集群、多核CPU或GPU硬件平台。

Claims (3)

1.一种加速克里金插值的方法，其特征在于它按以下步骤进行：

步骤一、在CPU端，按已知观测点的横坐标值对观测点从小到大排序；

步骤二、将排好序的数据传送到GPU，并启动GPU端克里金插值程序；

步骤三、在GPU中每个线程执行一个待估点的计算，所有线程首先搜索待估未知点x周围n个已知观测点x_j,j＝1,2,…,n，作为该未知点的邻近点，具体搜索过程如下：

1)将待估计点x横坐标值同已经排好序的观测点的横坐标比较，找到其在观测点中的排序，即x_i<x<x_i+1；

2)定义一大小为n的数组用于存储x的n个邻近点，以x点所在位置开始交替向左右两侧搜索，分两种情况处理：1、数组未被填满：将搜索到的点填入数组；2、数组已经被填满：如果搜索到的点满足下述条件1则替换数组内距x最大距离的邻近点，否则继续搜索；如果满足条件2或条件3则完成搜索；

其中，条件1为：x到当前搜索到的观测点的距离小于数组内所有存储的观测点到x的距离；

条件2为：x到当前搜索点的横坐标间的距离大于数组内所有存储的观测点到x的距离；

条件3为：所有观测点搜索完成；

步骤四、同步GPU中同一个workgroup内所有线程，待workgroup内所有的待估未知点的邻近点搜索完成后，执行下一步；

步骤五、所有workgroup内的第一个线程比较所有相邻待估未知点的邻近点，并标识具有相同邻近观测点的待估未知点；

其中，对于具有相同邻近点的待估点中最先出现的点，称为父节点；

步骤六、标记为父节点的GPU线程，根据变异函数公式及其拟合模型构造变异函数矩阵：

在这个步骤中，只有标记为父节点的线程参与了计算并保存结果矩阵，其它线程既没有消耗计算，也没有浪费存储空间；

步骤七、标记为父节点的GPU线程求K′的逆矩阵K′^-1；

步骤八、GPU内所有线程构造自己的矩阵：

步骤九、GPU内所有线程解方程组λ′＝M′·K′^-1，求出λ′，其中，μ是拉格朗日乘子；

步骤十、GPU内所有线程根据系数λ_i估计待测值x

2.如权利要求1所述的一种加速克里金插值的方法，其特征在于步骤五中具有相同邻近观测点的待估未知点的标识方法：每一个方格表示一个待估点，A、B…分别表示待估点的邻近点集合，如果相邻的待估点具有同样的邻近点，那么在方格中的字母也相同；

a、分别将每一个待估点同它的左侧的待估点比较，最左侧的待估点除外，如相同，在方格中表示为←；

b、分别将每一个待估点同它的上边的待估点相比较，最上层待估点除外，如相同，在方格中表示为↑；

c、对于具有相同邻近点的待估点中最先出现的点，从左向右，从上向下数，标为字母的待估点，称为父节点，并按从左到右，从上到下的顺序标记顺序号；

d、将其它非首次出现的待估点称为首次出现的待估点的子节点，将其填上其父节点的顺序号。

3.如权利要求2所述的一种加速克里金插值的方法，其特征在于该 方法适合于集群、多核CPU或GPU硬件平台。