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

王冰怡， 鄧維科， 余先川

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

?

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

王冰怡， 鄧維科， 余先川

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

條件模擬是一種計(jì)算非常耗時(shí)的高精度三維插值算法。針對(duì)串行條件模擬算法計(jì)算時(shí)間過(guò)長(zhǎng)的問(wèn)題，提出基于GPU的并行條件模擬算法，并進(jìn)行儲(chǔ)量估算。對(duì)條件模擬算法進(jìn)行并行分析，利用GPU的高度并行性，構(gòu)建CUDA通用計(jì)算開(kāi)發(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)域一門(mén)發(fā)展迅速的新興邊緣學(xué)科(王仁鐸等，1988)，在優(yōu)化采樣方案、處理不規(guī)則采樣及最優(yōu)化插值計(jì)算等方面具有明顯的優(yōu)點(diǎn)。條件模擬是地質(zhì)統(tǒng)計(jì)學(xué)中迅速發(fā)展的一門(mén)新技術(shù)，是一種新的蒙特卡羅(Monte Carlo)方法。采用條件模擬方法得出的模擬值，不僅能與已知數(shù)據(jù)Z(x)保持一致的數(shù)學(xué)期望、方差和分布函數(shù)，而且還能保持一致的協(xié)方差函數(shù)C(h)或變差函數(shù)γ(h)，同時(shí)還滿(mǎn)足在各已知點(diǎn)處的模擬值等于該點(diǎn)實(shí)測(cè)值的條件。條件模擬充分考慮了原始量的空間結(jié)構(gòu)性、隨機(jī)性和局部變異性，較好地解決了克立格法中難以克服的平滑效應(yīng)等問(wèn)題(李星等，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)是滿(mǎn)足二階平穩(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上滿(mǎn)足模擬值等于實(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í)間過(guò)長(zhǎng)，無(wú)法滿(mǎn)足實(shí)時(shí)分析的要求(白樹(shù)仁等，2011)，并行算法是解決這一問(wèn)題的一種有效方法(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ì)算開(kāi)發(fā)環(huán)境，提出了基于GPU的并行條件模擬算法，在不失真的情況下降低了運(yùn)算時(shí)間，提高了估算精確度。

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

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

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

統(tǒng)一計(jì)算設(shè)備架構(gòu)(Compute Unified Device Architecture，CUDA)是一個(gè)軟硬件協(xié)同的完整的解決方案。這種架構(gòu)可以使用GPU處理復(fù)雜的科學(xué)計(jì)算問(wèn)題，特別是極大數(shù)據(jù)量的并行計(jì)算問(wèn)題。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ǔ)量，其串行過(guò)程包括以下步驟。

(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ì)算，因此可以通過(guò)并行計(jì)算的思想來(lái)開(kāi)發(fā)。步驟(5)為加法運(yùn)算，可用CUDA的規(guī)約方式實(shí)現(xiàn)并行優(yōu)化。分析得到并行條件模擬插值算法步驟如下。

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

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

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

(4) GPU各線(xiàn)程：對(duì)于線(xiàn)程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ǔ)量，見(jiàn)公式(5)。

(5) GPU各線(xiàn)程：

① _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。過(guò)程(5)將計(jì)算累積結(jié)果的算法時(shí)間復(fù)雜度從O(n)降至O(logn)，并且通過(guò)GPU并行化計(jì)算，提高了運(yùn)算效率。

(6) CPU主線(xiàn)程：保存條件模擬結(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線(xiàn)程)；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ù)來(lái)源于中國(guó)地質(zhì)科學(xué)院礦產(chǎn)資源研究所

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

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

圖1　 CPU 與GPU 條件模擬耗時(shí)增長(zhǎng)曲線(xiàn)比較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的條件模擬過(guò)程耗時(shí)越長(zhǎng)，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)，因此通過(guò)表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的加速效果將越來(lái)越好，并行計(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)算過(guò)程。

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

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

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

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

王仁鐸，胡光道，1988. 線(xiàn)性地質(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；編輯：陳露

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

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

P612; P628+.5

A

1674-3636(2016)03-0507-05