王冰怡， 鄧維科， 余先川

(北京師范大學(xué)信息科學(xué)與技術(shù)學(xué)院，北京100875)

?

基于GPU的并行條件模擬算法及其在儲(chǔ)量估算中的應(yīng)用

王冰怡， 鄧維科， 余先川

(北京師范大學(xué)信息科學(xué)與技術(shù)學(xué)院，北京100875)

條件模擬是一種計(jì)算非常耗時(shí)的高精度三維插值算法。針對(duì)串行條件模擬算法計(jì)算時(shí)間過長的問題，提出基于GPU的并行條件模擬算法，并進(jìn)行儲(chǔ)量估算。對(duì)條件模擬算法進(jìn)行并行分析，利用GPU的高度并行性，構(gòu)建CUDA通用計(jì)算開發(fā)環(huán)境，實(shí)現(xiàn)串行條件模擬算法到并行條件模擬算法的轉(zhuǎn)換，使條件模擬算法的時(shí)間復(fù)雜度從O(n)降至O(logn)。并對(duì)西藏甲瑪銅礦進(jìn)行了儲(chǔ)量估算。實(shí)驗(yàn)結(jié)果表明，在安裝普通NVIDIA顯卡的計(jì)算機(jī)以及估算精度不下降的情況下，GPU并行條件模擬的計(jì)算效率比CPU串行條件模擬的計(jì)算效率提高了60倍以上。

條件模擬；GPU；并行計(jì)算；儲(chǔ)量估算；甲瑪銅礦；西藏

0　引　言

地質(zhì)統(tǒng)計(jì)學(xué)是數(shù)學(xué)地質(zhì)領(lǐng)域一門發(fā)展迅速的新興邊緣學(xué)科(王仁鐸等，1988)，在優(yōu)化采樣方案、處理不規(guī)則采樣及最優(yōu)化插值計(jì)算等方面具有明顯的優(yōu)點(diǎn)。條件模擬是地質(zhì)統(tǒng)計(jì)學(xué)中迅速發(fā)展的一門新技術(shù)，是一種新的蒙特卡羅(Monte Carlo)方法。采用條件模擬方法得出的模擬值，不僅能與已知數(shù)據(jù)Z(x)保持一致的數(shù)學(xué)期望、方差和分布函數(shù)，而且還能保持一致的協(xié)方差函數(shù)C(h)或變差函數(shù)γ(h)，同時(shí)還滿足在各已知點(diǎn)處的模擬值等于該點(diǎn)實(shí)測(cè)值的條件。條件模擬充分考慮了原始量的空間結(jié)構(gòu)性、隨機(jī)性和局部變異性，較好地解決了克立格法中難以克服的平滑效應(yīng)等問題(李星等，2003)。目前，條件模擬基本方法主要有轉(zhuǎn)向帶法(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)等。其中，誤差模擬方法作為基礎(chǔ)方法運(yùn)用最為廣泛，其主要思路為：在克立格估值的基礎(chǔ)上，人為增加一個(gè)被克立格法的平滑效應(yīng)濾掉的隨機(jī)誤差部分。

1　串行條件模擬算法及其應(yīng)用缺陷

設(shè)Z(x)是滿足二階平穩(wěn)假設(shè)的局域化變量，E[Z(x)]=m并存在協(xié)方差函數(shù)C(h)或變差函數(shù)γ(h)。要求Z(x)的條件模擬為ZCS(x)，即要求獲得與Z(x)同構(gòu)的區(qū)域化變量ZCS(x)的一個(gè)現(xiàn)實(shí)，且在已知點(diǎn)xa上滿足模擬值等于實(shí)測(cè)值。即：

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

(1)

(2)

(3)

(4)

式(4)為計(jì)算誤差條件模擬的實(shí)用型公式。根據(jù)公式(4)可知，首先求出非條件模擬值ZS(x)，再對(duì)實(shí)測(cè)點(diǎn)xa上的誤差值[Z(xa)-ZS(xa)],a=1,2,…,n進(jìn)行克立格估算，最后將二者相加得到ZCS(x)，此方法可以減少1次克立格方程組計(jì)算。

隨著已知信息的增加，估算層數(shù)眾多、網(wǎng)格劃分細(xì)的大型區(qū)域儲(chǔ)量時(shí)，串行條件模擬算法計(jì)算時(shí)間過長，無法滿足實(shí)時(shí)分析的要求(白樹仁等，2011)，并行算法是解決這一問題的一種有效方法(Pesquer et al., 2011；Hu et al., 2015；Liu et al., 2016)。利用GPU在處理通用計(jì)算時(shí)的高度并行性、優(yōu)秀的浮點(diǎn)運(yùn)算能力，構(gòu)建CUDA通用計(jì)算開發(fā)環(huán)境，提出了基于GPU的并行條件模擬算法，在不失真的情況下降低了運(yùn)算時(shí)間，提高了估算精確度。

2　基于GPU的并行條件模擬算法

2.1GPU與CUDA開發(fā)環(huán)境

GPU專注于解決數(shù)據(jù)并行計(jì)算的問題，具有極高的計(jì)算密度。GPU 并不要求對(duì)精密流的控制，可通過高速運(yùn)算來隱藏存儲(chǔ)器的訪問延遲，而不必使用較大的數(shù)據(jù)緩存。

統(tǒng)一計(jì)算設(shè)備架構(gòu)(Compute Unified Device Architecture，CUDA)是一個(gè)軟硬件協(xié)同的完整的解決方案。這種架構(gòu)可以使用GPU處理復(fù)雜的科學(xué)計(jì)算問題，特別是極大數(shù)據(jù)量的并行計(jì)算問題。CUDA 在GPU 架構(gòu)上將晶體管更多地投入到數(shù)據(jù)處理，減少數(shù)據(jù)緩存和流量控制對(duì)晶體管資源的消耗。

2.2并行條件模擬算法

并行條件模擬算法首先需要將地質(zhì)實(shí)體網(wǎng)格化，將其劃分為若干規(guī)則的網(wǎng)格塊段。條件模擬法依據(jù)空間已知點(diǎn)，對(duì)待模擬網(wǎng)格進(jìn)行條件模擬，并將所有待模擬網(wǎng)格相加得到礦體總儲(chǔ)量，其串行過程包括以下步驟。

(1) 輸入：空間已知點(diǎn)集合A，集合A元素為Pi(x,y,z,Z(x,y,z)),i=1,2,…,N；各向異性軸信息；輸出：條件模擬結(jié)果，即若干網(wǎng)格。

(2) 根據(jù)各向異性軸，對(duì)礦體實(shí)體進(jìn)行網(wǎng)格化，將地質(zhì)實(shí)體劃分為S個(gè)待模擬網(wǎng)格。

(3) fori= 1 toS

① 根據(jù)3個(gè)各向異性軸信息，形成1個(gè)橢球形的搜索區(qū)域；

② 以網(wǎng)格Gridi為橢球中心，搜索落在此橢球搜索區(qū)域內(nèi)、距離該中心最近的若干空間已知點(diǎn)，將其保存為集合B；

③ 根據(jù)B中保存的已知點(diǎn)，計(jì)算待模擬網(wǎng)格Gridi的條件模擬結(jié)果ZSi(x,y,z)；

④ 若條件模擬為序貫方法(序貫高斯條件模擬與序貫指示條件模擬)，則需要將條件模擬結(jié)果放入已知點(diǎn)集合中；

end for。

(4) 根據(jù)每個(gè)網(wǎng)格儲(chǔ)量估算公式計(jì)算儲(chǔ)量結(jié)果Wi：

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

(5)

式(5)中，Wi為儲(chǔ)量估算結(jié)果，x,y,z為網(wǎng)格塊段位置，ZSi(x,y,z)為該網(wǎng)格處的條件模擬值，Vi(x,y,z)為網(wǎng)格塊段體積，p為礦石體重(即密度)。

(5) 將所有待模擬網(wǎng)格計(jì)算后的儲(chǔ)量，相加得到礦體總條件模擬量W：

(6)

對(duì)于上述步驟(3)，每個(gè)待模擬的網(wǎng)格都是一個(gè)獨(dú)立的計(jì)算，因此可以通過并行計(jì)算的思想來開發(fā)。步驟(5)為加法運(yùn)算，可用CUDA的規(guī)約方式實(shí)現(xiàn)并行優(yōu)化。分析得到并行條件模擬插值算法步驟如下。

(1) CPU主線程：同串行條件模擬算法中的步驟(1)。

(2) CPU主線程：執(zhí)行串行條件模擬算法中的步驟(2)，并將網(wǎng)格信息作為參數(shù)輸入CUDA的kernel核函數(shù)。

(3) CPU主線程：設(shè)置線程屬性，包括CUDA計(jì)算環(huán)境的Grid、Block中線程數(shù)。

(4) GPU各線程：對(duì)于線程i，計(jì)算網(wǎng)格xi對(duì)應(yīng)的插值結(jié)果Zi(x,y,z)。

① 利用變異性分析結(jié)果，獲取3個(gè)各向異性軸信息，形成一個(gè)橢球形的搜索區(qū)域。

② 以網(wǎng)格i為橢球中心，搜索落在此橢球搜索區(qū)域內(nèi)的空間已知點(diǎn)，將其保存。

③ 根據(jù)步驟②中保存的已知點(diǎn)，利用條件模擬法進(jìn)行計(jì)算得到待模擬網(wǎng)格i的條件模擬值ZSi(x,y,z)。

④ 根據(jù)每個(gè)網(wǎng)格條件模擬后的結(jié)果，計(jì)算儲(chǔ)量，見公式(5)。

(5) GPU各線程：

① _syncthreads()； //CUDA同步函數(shù)

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

③ while (i != 1)

if (j < i)

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

end if；

_syncthreads()；//CUDA同步函數(shù)

i = i/ 2；

end while；

④ return W[1]；

其中，W[1]為相加后的礦體條件模擬結(jié)果W。過程(5)將計(jì)算累積結(jié)果的算法時(shí)間復(fù)雜度從O(n)降至O(logn)，并且通過GPU并行化計(jì)算，提高了運(yùn)算效率。

(6) CPU主線程：保存條件模擬結(jié)果W。

3　實(shí)驗(yàn)設(shè)置與結(jié)果分析

3.1實(shí)驗(yàn)設(shè)置

為了驗(yàn)證基于GPU的并行條件模擬算法的效率，在Windows操作系統(tǒng)下，使用GPU方法與傳統(tǒng)CPU方法分別進(jìn)行條件模擬實(shí)驗(yàn)，其計(jì)算環(huán)境如下。

GPU運(yùn)算平臺(tái)為NVIDIA GeForce GT445M顯卡(144 CUDA核)；GPU最低主頻1.18 GHz，動(dòng)態(tài)調(diào)頻；顯存3 072 M(DDR3 800 MHz)，位寬192 bit；CUDA驅(qū)動(dòng)版本4.10/4.0；CUDA運(yùn)算版本2.1。操作系統(tǒng)Windows 7(64 bit)。

CPU運(yùn)算平臺(tái)為Intel i7 Q740(4核8線程)；CPU最低主頻1.73 GHz，動(dòng)態(tài)調(diào)頻；內(nèi)存8 G(DDR3 1 333 MHz)。操作系統(tǒng)Windows 7(64 bit)。

采用西藏甲瑪銅礦實(shí)驗(yàn)數(shù)據(jù)，鉆孔采樣點(diǎn)統(tǒng)計(jì)信息如表1所示。

注：數(shù)據(jù)來源于中國地質(zhì)科學(xué)院礦產(chǎn)資源研究所

3.2實(shí)驗(yàn)結(jié)果分析

條件模擬算法在CPU平臺(tái)和GPU平臺(tái)上的計(jì)算時(shí)間增長曲線如圖1所示。從圖中看出，隨著條件模擬網(wǎng)格數(shù)的增加，條件模擬的CPU計(jì)算時(shí)間越來越長，GPU的效率依然很高，優(yōu)越性越明顯。

圖1　 CPU 與GPU 條件模擬耗時(shí)增長曲線比較Fig.1　Comparison of time-consuming growth curves of CPU and GPU conditional simulation

條件模擬在CPU 平臺(tái)與GPU 平臺(tái)的運(yùn)算時(shí)間比較如表2所示。隨著網(wǎng)格精細(xì)程度的增加，基于CPU的條件模擬過程耗時(shí)越長，GPU的加速效果越明顯，并在達(dá)到峰值后保持平穩(wěn)。

表2　CPU與GPU 的條件模擬運(yùn)算時(shí)間比較

分別利用基于CPU的條件模擬和基于GPU的條件模擬，對(duì)西藏甲瑪銅礦進(jìn)行儲(chǔ)量估算(表3)。理論上，由于條件模擬具有隨機(jī)性，GPU平臺(tái)的條件模擬結(jié)果與CPU平臺(tái)的條件模擬結(jié)果會(huì)有所差別，但均在克立格方差的可控范圍內(nèi)，因此通過表3的對(duì)比分析，也同樣驗(yàn)證了基于GPU將條件模擬并行化的準(zhǔn)確性與穩(wěn)定性。

表3　CPU與GPU條件模擬的儲(chǔ)量估算結(jié)果比較

分析圖2的GPU加速比，可以看出：當(dāng)數(shù)據(jù)數(shù)量較少時(shí)，GPU的加速效果不明顯；隨著待模擬網(wǎng)格數(shù)量的增加，計(jì)算時(shí)間也隨之增加，GPU的加速效果將越來越好，并行計(jì)算的優(yōu)越性越明顯；直到加速至GPU的并行加速極限時(shí)，加速比趨于平穩(wěn)。

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

4　結(jié)　論

(1) 在已有的條件模擬模型基礎(chǔ)上，針對(duì)三維實(shí)體模型的復(fù)雜性，提出分布式算法，改進(jìn)了運(yùn)算過程。

(2) 提出了基于GPU并行計(jì)算的并行條件模擬算法，將時(shí)間復(fù)雜度從O(n)降至O(logn)，并進(jìn)行了儲(chǔ)量估算。

(3) 實(shí)驗(yàn)結(jié)果表明，提出的并行條件模擬在估計(jì)精度不變或不下降的情況下具有很高的加速比與并行計(jì)算效率，適用于固體礦產(chǎn)儲(chǔ)量計(jì)算。

白樹仁，李濤，寧錦陽，2011. 基于MPI的并行Kriging空間降水插值[J]. 計(jì)算技術(shù)與自動(dòng)化，30(1)：71-74.

李星，趙彥超，王國慶,2003. 條件模擬的原理及應(yīng)用[J]. 地球?qū)W報(bào), 24(增刊1)：226-228.

王仁鐸，胡光道，1988. 線性地質(zhì)統(tǒng)計(jì)學(xué)[M]. 北京：地質(zhì)出版社.

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；編輯：陳露

國土資源部公益性行業(yè)科研專項(xiàng)經(jīng)費(fèi)項(xiàng)目(201511079-02)，國家自然科學(xué)基金項(xiàng)目(41272359)，教育部博士點(diǎn)基金項(xiàng)目(20120003110032)

王冰怡(1995—)，女，碩士研究生，計(jì)算機(jī)軟件與理論專業(yè)，主要研究方向?yàn)榭臻g信息數(shù)據(jù)處理與挖掘，E-mail：201521210034@mail.bnu.edu.com

P612; P628+.5

A

1674-3636(2016)03-0507-05