基于GPU的并行条件模拟算法及其在储量估算中的应用

王冰怡， 邓维科， 余先川

(北京师范大学信息科学与技术学院，北京100875)



基于GPU的并行条件模拟算法及其在储量估算中的应用

王冰怡， 邓维科， 余先川

(北京师范大学信息科学与技术学院，北京100875)

条件模拟是一种计算非常耗时的高精度三维插值算法。针对串行条件模拟算法计算时间过长的问题，提出基于GPU的并行条件模拟算法，并进行储量估算。对条件模拟算法进行并行分析，利用GPU的高度并行性，构建CUDA通用计算开发环境，实现串行条件模拟算法到并行条件模拟算法的转换，使条件模拟算法的时间复杂度从O(n)降至O(logn)。并对西藏甲玛铜矿进行了储量估算。实验结果表明，在安装普通NVIDIA显卡的计算机以及估算精度不下降的情况下，GPU并行条件模拟的计算效率比CPU串行条件模拟的计算效率提高了60倍以上。

条件模拟；GPU；并行计算；储量估算；甲玛铜矿；西藏

0　引　言

地质统计学是数学地质领域一门发展迅速的新兴边缘学科(王仁铎等，1988)，在优化采样方案、处理不规则采样及最优化插值计算等方面具有明显的优点。条件模拟是地质统计学中迅速发展的一门新技术，是一种新的蒙特卡罗(Monte Carlo)方法。采用条件模拟方法得出的模拟值，不仅能与已知数据Z(x)保持一致的数学期望、方差和分布函数，而且还能保持一致的协方差函数C(h)或变差函数γ(h)，同时还满足在各已知点处的模拟值等于该点实测值的条件。条件模拟充分考虑了原始量的空间结构性、随机性和局部变异性，较好地解决了克立格法中难以克服的平滑效应等问题(李星等，2003)。目前，条件模拟基本方法主要有转向带法(Brooker，1985；Dean，1995)、误差模拟方法(Boudreault et al., 2016；Liang et al., 2016)、傅立叶变换法(Pinheiro et al., 2016)和序贯条件模拟方法(Juang et al., 2004；Sojdehee et al., 2015；Yamamoto et al., 2015；Ribeiro et al., 2016)等。其中，误差模拟方法作为基础方法运用最为广泛，其主要思路为：在克立格估值的基础上，人为增加一个被克立格法的平滑效应滤掉的随机误差部分。

1　串行条件模拟算法及其应用缺陷

设Z(x)是满足二阶平稳假设的局域化变量，E[Z(x)]=m并存在协方差函数C(h)或变差函数γ(h)。要求Z(x)的条件模拟为ZCS(x)，即要求获得与Z(x)同构的区域化变量ZCS(x)的一个现实，且在已知点xa上满足模拟值等于实测值。即：

ZCS(xa)=Z(xa),∀xaa=1,2,…,n

(1)

(2)

(3)

(4)

式(4)为计算误差条件模拟的实用型公式。根据公式(4)可知，首先求出非条件模拟值ZS(x)，再对实测点xa上的误差值[Z(xa)-ZS(xa)],a=1,2,…,n进行克立格估算，最后将二者相加得到ZCS(x)，此方法可以减少1次克立格方程组计算。

随着已知信息的增加，估算层数众多、网格划分细的大型区域储量时，串行条件模拟算法计算时间过长，无法满足实时分析的要求(白树仁等，2011)，并行算法是解决这一问题的一种有效方法(Pesquer et al., 2011；Hu et al., 2015；Liu et al., 2016)。利用GPU在处理通用计算时的高度并行性、优秀的浮点运算能力，构建CUDA通用计算开发环境，提出了基于GPU的并行条件模拟算法，在不失真的情况下降低了运算时间，提高了估算精确度。

2　基于GPU的并行条件模拟算法

2.1GPU与CUDA开发环境

GPU专注于解决数据并行计算的问题，具有极高的计算密度。GPU 并不要求对精密流的控制，可通过高速运算来隐藏存储器的访问延迟，而不必使用较大的数据缓存。

统一计算设备架构(Compute Unified Device Architecture，CUDA)是一个软硬件协同的完整的解决方案。这种架构可以使用GPU处理复杂的科学计算问题，特别是极大数据量的并行计算问题。CUDA 在GPU 架构上将晶体管更多地投入到数据处理，减少数据缓存和流量控制对晶体管资源的消耗。

2.2并行条件模拟算法

并行条件模拟算法首先需要将地质实体网格化，将其划分为若干规则的网格块段。条件模拟法依据空间已知点，对待模拟网格进行条件模拟，并将所有待模拟网格相加得到矿体总储量，其串行过程包括以下步骤。

(1) 输入：空间已知点集合A，集合A元素为Pi(x,y,z,Z(x,y,z)),i=1,2,…,N；各向异性轴信息；输出：条件模拟结果，即若干网格。

(2) 根据各向异性轴，对矿体实体进行网格化，将地质实体划分为S个待模拟网格。

(3) fori= 1 toS

① 根据3个各向异性轴信息，形成1个椭球形的搜索区域；

② 以网格Gridi为椭球中心，搜索落在此椭球搜索区域内、距离该中心最近的若干空间已知点，将其保存为集合B；

③ 根据B中保存的已知点，计算待模拟网格Gridi的条件模拟结果ZSi(x,y,z)；

④ 若条件模拟为序贯方法(序贯高斯条件模拟与序贯指示条件模拟)，则需要将条件模拟结果放入已知点集合中；

end for。

(4) 根据每个网格储量估算公式计算储量结果Wi：

Wi=ZSi(x,y,z)·p·Vi(x,y,z)/100

(5)

式(5)中，Wi为储量估算结果，x,y,z为网格块段位置，ZSi(x,y,z)为该网格处的条件模拟值，Vi(x,y,z)为网格块段体积，p为矿石体重(即密度)。

(5) 将所有待模拟网格计算后的储量，相加得到矿体总条件模拟量W：

(6)

对于上述步骤(3)，每个待模拟的网格都是一个独立的计算，因此可以通过并行计算的思想来开发。步骤(5)为加法运算，可用CUDA的规约方式实现并行优化。分析得到并行条件模拟插值算法步骤如下。

(1) CPU主线程：同串行条件模拟算法中的步骤(1)。

(2) CPU主线程：执行串行条件模拟算法中的步骤(2)，并将网格信息作为参数输入CUDA的kernel核函数。

(3) CPU主线程：设置线程属性，包括CUDA计算环境的Grid、Block中线程数。

(4) GPU各线程：对于线程i，计算网格xi对应的插值结果Zi(x,y,z)。

① 利用变异性分析结果，获取3个各向异性轴信息，形成一个椭球形的搜索区域。

② 以网格i为椭球中心，搜索落在此椭球搜索区域内的空间已知点，将其保存。

③ 根据步骤②中保存的已知点，利用条件模拟法进行计算得到待模拟网格i的条件模拟值ZSi(x,y,z)。

④ 根据每个网格条件模拟后的结果，计算储量，见公式(5)。

(5) GPU各线程：

① _syncthreads()； //CUDA同步函数

② i= blockDim.x/ 2；j = 1；

③ while (i != 1)

if (j < i)

then W[i] = W[i] + W[j + i]；

end if；

_syncthreads()；//CUDA同步函数

i = i/ 2；

end while；

④ return W[1]；

其中，W[1]为相加后的矿体条件模拟结果W。过程(5)将计算累积结果的算法时间复杂度从O(n)降至O(logn)，并且通过GPU并行化计算，提高了运算效率。

(6) CPU主线程：保存条件模拟结果W。

3　实验设置与结果分析

3.1实验设置

为了验证基于GPU的并行条件模拟算法的效率，在Windows操作系统下，使用GPU方法与传统CPU方法分别进行条件模拟实验，其计算环境如下。

GPU运算平台为NVIDIA GeForce GT445M显卡(144 CUDA核)；GPU最低主频1.18 GHz，动态调频；显存3 072 M(DDR3 800 MHz)，位宽192 bit；CUDA驱动版本4.10/4.0；CUDA运算版本2.1。操作系统Windows 7(64 bit)。

CPU运算平台为Intel i7 Q740(4核8线程)；CPU最低主频1.73 GHz，动态调频；内存8 G(DDR3 1 333 MHz)。操作系统Windows 7(64 bit)。

采用西藏甲玛铜矿实验数据，钻孔采样点统计信息如表1所示。

注：数据来源于中国地质科学院矿产资源研究所

3.2实验结果分析

条件模拟算法在CPU平台和GPU平台上的计算时间增长曲线如图1所示。从图中看出，随着条件模拟网格数的增加，条件模拟的CPU计算时间越来越长，GPU的效率依然很高，优越性越明显。

图1　 CPU 与GPU 条件模拟耗时增长曲线比较Fig.1　Comparison of time-consuming growth curves of CPU and GPU conditional simulation

条件模拟在CPU 平台与GPU 平台的运算时间比较如表2所示。随着网格精细程度的增加，基于CPU的条件模拟过程耗时越长，GPU的加速效果越明显，并在达到峰值后保持平稳。

表2　CPU与GPU 的条件模拟运算时间比较

分别利用基于CPU的条件模拟和基于GPU的条件模拟，对西藏甲玛铜矿进行储量估算(表3)。理论上，由于条件模拟具有随机性，GPU平台的条件模拟结果与CPU平台的条件模拟结果会有所差别，但均在克立格方差的可控范围内，因此通过表3的对比分析，也同样验证了基于GPU将条件模拟并行化的准确性与稳定性。

表3　CPU与GPU条件模拟的储量估算结果比较

分析图2的GPU加速比，可以看出：当数据数量较少时，GPU的加速效果不明显；随着待模拟网格数量的增加，计算时间也随之增加，GPU的加速效果将越来越好，并行计算的优越性越明显；直到加速至GPU的并行加速极限时，加速比趋于平稳。

图2　 GPU的加速比分析Fig.2　Analysis of speed-up ratio for GPU

4　结　论

(1) 在已有的条件模拟模型基础上，针对三维实体模型的复杂性，提出分布式算法，改进了运算过程。

(2) 提出了基于GPU并行计算的并行条件模拟算法，将时间复杂度从O(n)降至O(logn)，并进行了储量估算。

(3) 实验结果表明，提出的并行条件模拟在估计精度不变或不下降的情况下具有很高的加速比与并行计算效率，适用于固体矿产储量计算。

白树仁，李涛，宁锦阳，2011. 基于MPI的并行Kriging空间降水插值[J]. 计算技术与自动化，30(1)：71-74.

李星，赵彦超，王国庆,2003. 条件模拟的原理及应用[J]. 地球学报, 24(增刊1)：226-228.

王仁铎，胡光道，1988. 线性地质统计学[M]. 北京：地质出版社.

BROOKER P I, 1985. Two-dimensional simulation by turning bands[J]. Journal of the International Association for Mathematical Geology, 17(1): 81-90.

BOUDREAULT J P, DUBÉ J S, MARCOTTE D, 2016. Quantification and minimization of uncertainty by geostatistical simulations during the characterization of contaminated sites: 3-D approach to a multi-element contamination[J]. Geoderma, 264: 214-226.

DEAN S O,1995. Moving averages for Gaussian simulation in two and three dimensions[J]. Mathematical Geology, 27(8): 939-960.

HU H D, SHU H, 2015. An improved coarse-grained parallel algorithm for computational acceleration of ordinary Kriging interpolation[J]. Computers & Geosciences, 78(2): 44-52.JUANG K W, CHEN Y S, LEE D Y, 2004. Using sequential indicator simulation to assess the uncertainty of delineating heavy-metal contaminated soils[J]. Environmental Pollution, 127(2): 229-238.

LIANG M, MARCOTTE D, SHAMSIPOUR P, 2016. Simulation of non-linear coregionalization models by FFTMA[J]. Computers & Geosciences, 89(1): 220-231.

LIU L, WU C L, WANG Z B, 2016. Parallelization of the Kriging Algorithm in stochastic simulation with GPU accelerators[C]//BIAN F, XIE Y C. Geo-Informatics in Resource Management and Sustainable Ecosystem. Berlin: Springer, 197-205.

OLIVER D S, 1995. Moving averages for Gaussian simulation in two and three dimensions[J]. Mathematical Geology, 27(8): 939-960.

PESQUER L, CORTÉS A, PONS X, 2011. Parallel ordinary Kriging interpolation incorporating automatic variogram fitting[J]. Computers & Geosciences, 37(4): 464-473.

PINHEIRO M, VALLEJOS J, MIRANDA T, et al., 2016. Geostatistical simulation to map the spatial heterogeneity of geomechanical parameters: A case study with rock mass rating[J]. Engineering Geology, 205: 93-103.

RIBEIRO S, CAINETA J COSTA A C, et al., 2016. Detection of inhomogeneities in precipitation time series in Portugal using direct sequential simulation[J]. Atmospheric Research, 171: 147-158.

SOJDEHEE M, RASA I, NEZAFATI N, et al., 2015. Probabilistic modeling of mineralized zones in Daralu copper deposit (SE Iran) using sequential indicator simulation[J]. Arabian Journal of Geosciences, 8(10): 8449-8459.

YAMAMOTO J K, LANDIM P M B, KIKUDA A T, et al., 2015. Post-processing of sequential indicator simulation realizations for modeling geologic bodies[J]. Computational Geosciences, 19(1): 257-266.

The parallel conditional simulation algorithm based on GPU and its application in reserve estimation

WANG Bingyi, DENG Weike, YU Xianchuan

(College of Information Science and Technology, Beijing Normal University, Beijing 100875, China)

Conditional simulation is a very time-consuming algorithm of three-dimensional high-precision interpolation. To solve the time-consuming problem in serial conditional simulation algorithm, we proposed a parallel conditional simulation algorithm based on GPU, and applied it to ore reserve estimation. First, this work conducted parallel analysis on the serial conditional simulation algorithm, and constructed CUDA general computing development environment based on the high parallelism of GPU. This realized the transformation of serial conditional simulation algorithm to parallel conditional simulation algorithm, and cut the time complexity of conditional simulation algorithm fromO(n) toO(logn). Experimental results of the Jiama copper deposit in Tibet on a computer with NVIDIA graphics demonstrated that, compared with the serial conditional algorithm, the GPU parallel conditional simulation algorithm could improve the efficiency by more than 60 times, without estimation accuracy loss.

conditional simulation; GPU; parallel computing; reserve estimation; Jiama copper deposit; Tibet

10.3969/j.issn.1674-3636.2016.03.507

2016-06-02；

2016-06-22；编辑：陈露

国土资源部公益性行业科研专项经费项目(201511079-02)，国家自然科学基金项目(41272359)，教育部博士点基金项目(20120003110032)

王冰怡(1995—)，女，硕士研究生，计算机软件与理论专业，主要研究方向为空间信息数据处理与挖掘，E-mail：201521210034@mail.bnu.edu.com

P612; P628+.5

A

1674-3636(2016)03-0507-05