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

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

【技术实现步骤摘要】

本专利技术涉及地学领域，尤其涉及一种加速克里金插值的方法。
技术介绍
克里金插值算法也称为空间自协方差最佳插值法，它最早由南非矿业工程师D.G.Krige提出，因此也以他的名字命名该算法。克里金插值算法广泛用于地学、环境科学等领域。它是根据未知点周围一定范围内观测点的值来估计未知点取值的一种方法。该方法是一种线性无偏方法，内插结果准确率高。然而在实际应用中，待插值地区面往往过大，插值密度也要求很高，这就造成了克里金插值速度很慢，严重限制了它的应用。随着并行计算技术的发展，基于这些新技术加速克里金算法成为一个研究热点。目前，克里金加速算法主要分为以下几个方面：一、基于集群技术集群由大量计算机通过高速网络连接组成，普通的商用集群一般包含上百个节点。集群的特点是相对普通PC(个人电脑)，运算能力比较强，但也存在着能耗高，管理成本高，应用难度高的缺点。基于集群硬件平台来加速克里金插值方法通常使用消息传递接口(MPI)，具体实现方法包括两种：一是将算法中每一个可以并行的部分在集群中以并行方式完成；二是以每个未知点的估计做为一个并行单元在集群上并行执行。前者编程相对容易，因为克里金插值中的可并行部分是矩阵计算，现有大量专业人员编写的可在集群上并行运行的代码库可以使用，但是这种并行方式对集群的利用率不高，同时会产生大量的节点间通讯，对速度有很大的影响；后者能充分利用集群的计算资源，但是需要编写克里金算法中的包括矩阵计算的各种数值算法，因此实现成本高，对编写代码人员的编程能力要求较高。二、基于多核CPU由于工艺原因，CPU的核心频率遇到了瓶颈，CPU的主流发展朝多核方向发展，充分利用多核CPU的资源也成为加速克里金的一种主要方法。多核CPU成本较集群低，实现相对容易，节点间通讯速度也较快，但受CPU核心数量的限制，其能够并行执行的线程数有限，因此加速比潜力不如集群或GPU。基于多核CPU加速克里金插值一般使用OpenMP框架，具体的加速方法也分为两类：一是将可并行执行的矩阵计算在多核CPU上并行执行；二是以每个待估的未知点做为基本的并行单元完成克里金插值的并行执行。三、基于GPUGPU平台是一种有前景的并行加速平台。GPU相对集群成本低很多，能耗也小很多，而且GPU包含的可并行执行的计算核心数量远高于多核CPU，主流的GPU都拥有几千个计算核心。基于GPU硬件平台的加速技术主要基于两种软件平台：统一计算设备架构(CUDA)和开放运算语言(OpenCL)。具有加速方法同集群和多核CPU一样，分为两类：一是将可并行执行的矩阵计算在多核CPU上并行执行；二是以每个待估的未知点做为基本的并行单元完成克里金插值的并行执行。四、组合方法由于一种计算设备可能同时包含上述几种计算平台，因此上述几种加速算法也可以混合应用，以最大程度挖掘计算平台的潜力。
技术实现思路
克里金插值是地学领域经常用到的一种插值方法，同其它插值方法相比，例如：反距离加权插值法、最近邻点插值法、线性插值三角网法、局部多项式法等，其插值精度较高。但是克里金插值是一种严重消耗计算和内存资源的算法，尤其是在大面积区域中做高密度克金金插值是相当慢的，严重地影响了克里金插值的实际应用价值。本专利技术针对克金格插值算法提出一种新的解决方案，能极大的提高克里金插值的速度。本专利技术为解决现有加速克里金插值的方法仅从计算平台的计算资源的角度来加速，克里金插值中仍存在大量的冗余计算和存储的问题，而提出一种加速克里金插值的方法。一种加速克里金插值的方法，按以下步骤进行：步骤一、在CPU端，按已知观测点的横坐标值对观测点从小到大排序；步骤二、将排好序的数据传送到GPU，并启动GPU端克里金插值程序；步骤三、在GPU中每个线程执行一个待估点的计算，所有线程首先搜索待估未知点x周围n个已知观测点xj(j＝1,2,...,n)，作为该未知点的邻近点；步骤四、同步GPU中同一个workgroup内所有线程，待workgroup内所有的待估未知点的邻近点搜索完成后，执行下一步；步骤五、所有workgroup内的第一个线程比较所有相邻待估未知点的邻近点，并标识具有相同邻近观测点的待估未知点；其中，对于具有相同邻近点的待估点中最先出现的点，称为父节点；步骤六、标记为父节点的GPU线程，根据变异函数公式及其拟合模型构造变异函数矩阵：K′=γ(x1,x1)γ(x1,x2)...γ(x1,xn)1γ(x2,x1)γ(x2,x2)...γ(x2,xn)1...............γ(xn,x1)γ(xn,x2)...γ(xn,xn)111...10;]]>在这个步骤中，只有标记为父节点的线程参与了计算并保存结果矩阵，其它线程即没有消耗计算，也没有浪费存储空间；步骤七、标记为父节点的GPU线程求K′的逆矩阵K′-1；步骤八、GPU内所有线程构造自己的矩阵：M′=γ‾(x1,x)γ‾(x2,x)...γ‾(xn,x)1;]]>步骤九、GPU内所有线程解方程序组λ′＝M′×K′-1，求出λ′，其中，λ′=λ1λ2...λnμ,]]>μ是拉格朗日乘子；步骤十、GPU内所有线程根据系数λi估计待测值xZok*(x)=Σi=1nλiZ(xi).]]>本专利技术包括以下有益效果：1、传统方法仅从计算平台的计算资源的角度来加速，本专利技术方法从根本上解决了克里金插值法的运算速度慢的问题，在克里金插值过程中，选择离待估计最近的一定数量观测点来计算权值系数，相邻的待估点通常有相同的邻近已知观测点，由于待估未知点的取值完全取决于其邻近已知观测点的值，因此，对于具有相同邻近观测点的待估点，只需计算其中一个待估点，由此减少大量的冗余计算和中间计算过程中的矩阵存储量。2、本专利技术方法提出一种搜索邻近点的方法，减少了搜索范围，从而提高了搜索速度。附图说明图1为具有相同邻近观测点的待估未知点的标识过程示意图。具体实施方式为使本专利技术的上述目的、特征和优点能够更加明显易懂，下面结合图1和具体实施方式对本专利技术作进一步详细的说明。具体实施方式一、本实施方式所述的一种加速克里金插值的方法，按以下步骤进行：步骤一、在CPU端，按已知观测点本文档来自技高网...

【技术保护点】
一种加速克里金插值的方法，其特征在于它按以下步骤进行：步骤一、在CPU端，按已知观测点的横坐标值对观测点从小到大排序；步骤二、将排好序的数据传送到GPU，并启动GPU端克里金插值程序；步骤三、在GPU中每个线程执行一个待估点的计算，所有线程首先搜索待估未知点x周围n个已知观测点xj(j＝1,2,...,n)，作为该未知点的邻近点；步骤四、同步GPU中同一个workgroup内所有线程，待workgroup内所有的待估未知点的邻近点搜索完成后，执行下一步；步骤五、所有workgroup内的第一个线程比较所有相邻待估未知点的邻近点，并标识具有相同邻近观测点的待估未知点；其中，对于具有相同邻近点的待估点中最先出现的点，称为父节点；步骤六、标记为父节点的GPU线程，根据变异函数公式及其拟合模型构造变异函数矩阵：K′=γ(x1,x1)γ(x1,x2)...γ(x1,xn)1γ(x2,x1)γ(x2,x2)...γ(x2,xn)1...............γ(xn,x1)γ(xn,x2)...γ(xn,xn)111...10;]]>在这个步骤中，只有标记为父节点的线程参与了计算并保存结果矩阵，其它线程即没有消耗计算，也没有浪费存储空间；步骤七、标记为父节点的GPU线程求K′的逆矩阵K′‑1；步骤八、GPU内所有线程构造自己的矩阵：M′=γ‾(x1,x)γ‾(x2,x)...γ‾(xn,x)1;]]>步骤九、GPU内所有线程解方程组λ′＝M′×K′‑1，求出λ′，其中，λ′=λ1λ2...λnμ,]]>μ是拉格朗日乘子；步骤十、GPU内所有线程根据系数λi估计待测值xZok*(x)=Σi=1nλiZ(xi).]]>...

【技术特征摘要】
1.一种加速克里金插值的方法，其特征在于它按以下步骤进行：
步骤一、在CPU端，按已知观测点的横坐标值对观测点从小到大排序；
步骤二、将排好序的数据传送到GPU，并启动GPU端克里金插值程序；
步骤三、在GPU中每个线程执行一个待估点的计算，所有线程首先搜索待估未知点x周
围n个已知观测点xj(j＝1,2,...,n)，作为该未知点的邻近点；
步骤四、同步GPU中同一个workgroup内所有线程，待workgroup内所有的待估未知点的
邻近点搜索完成后，执行下一步；
步骤五、所有workgroup内的第一个线程比较所有相邻待估未知点的邻近点，并标识具
有相同邻近观测点的待估未知点；
其中，对于具有相同邻近点的待估点中最先出现的点，称为父节点；
步骤六、标记为父节点的GPU线程，根据变异函数公式及其拟合模型构造变异函数矩阵：K′=γ(x1,x1)γ(x1,x2)...γ(x1,xn)1γ(x2,x1)γ(x2,x2)...γ(x2,xn)1...............γ(xn,x1)γ(xn,x2)...γ(xn,xn)111...10;]]>在这个步骤中，只有标记为父节点的线程参与了计算并保存结果矩阵，其它线程即没
有消耗计算，也没有浪费存储空间；
步骤七、标记为父节点的GPU线程求K′的逆矩阵K′-1；
步骤八、GPU内所有线程构造自己的矩阵：M′=γ‾(x1,x)γ‾(x2,x)...γ‾(xn,x)1;]]>步骤九、GPU内所有线程解方程组λ′＝M′×K′-1，求出λ′，其中，λ′=λ1λ2...λnμ,]]>μ是拉格朗日
乘子；<...

【专利技术属性】
技术研发人员：姜春雷，张淑清，张策，
申请(专利权)人：中国科学院东北地理与农业生态研究所，
类型：发明
国别省市：黑龙江;23