羅德新，謝凱，夏巍，廖虎

（長江大學電子信息學院，湖北 荊州434023）

油氣資源是關(guān)系我國現(xiàn)代化建設(shè)全局和國家安全的重要戰(zhàn)略資源，在國民經(jīng)濟中也占有非常重要的地位。隨著經(jīng)濟的快速發(fā)展，世界各國都不同程度地面臨著油氣需求量急劇增加的問題，因此地質(zhì)勘探技術(shù)的發(fā)展，對準確高效的尋找油氣儲層起著至關(guān)重要的作用。三維地質(zhì)建模是指運用計算機技術(shù)，在三維環(huán)境下，將空間信息管理、地質(zhì)解釋、空間分析和預測、地學統(tǒng)計、實體內(nèi)容分析以及圖形可視化等工具結(jié)合起來以用于地質(zhì)研究。在復雜地質(zhì)體的環(huán)境下，地質(zhì)構(gòu)造復雜，巖性巖相不穩(wěn)定，地層信息難以獲取，這些問題都給地質(zhì)建模帶來很大的影響［1，2］。近年來，我國在油氣勘探方面的研究蓬勃發(fā)展，開發(fā)了具有自主知識產(chǎn)權(quán)的地質(zhì)建模系統(tǒng)［3，4］，其算法構(gòu)架和功能穩(wěn)定方面不乏先進性，但是在建模質(zhì)量、速度及準確度方面與國外同類型軟件相比存在明顯不足。

隨著對復雜地質(zhì)體的勘探不斷深入，復雜幾何模型的建模要求也越來越高。在描述極其復雜的地質(zhì)體時，經(jīng)常需要對不規(guī)則分布的大數(shù)據(jù)進行地質(zhì)建模，這些地質(zhì)模型經(jīng)常包含斷層、褶曲、開挖邊界等，其網(wǎng)格通常由成千上萬個多邊形面片組成.因為這種網(wǎng)格點分布具有很大的隨機性，多邊形構(gòu)造復雜，將多邊形映射為圖像的過程需要大量的計算，因而所需的計算時間相當長，數(shù)據(jù)存儲量也特別大，這對建模算法和建模速度提出了新的要求［5，6］。為此，筆者提出了一種基于GPU并行加速的快速建模方法。

1 建立規(guī)則的初始矩形網(wǎng)格模型

對于復雜地質(zhì)體的網(wǎng)格模型，首先在數(shù)據(jù)中查找4個邊界分布點，根據(jù)4個邊界分布點確定網(wǎng)格模型的邊界值。為了提高其建模的精確度，通過數(shù)據(jù)中采樣點的平均間距來確定網(wǎng)格的劃分間距。采樣點的X方向的平均間距為x，Y方向的平均間距為y，網(wǎng)格X方向劃分間距為Δx，Y方向劃分間距為Δy，劃分網(wǎng)格如圖1所示，則每個網(wǎng)格單元中平均含有采樣點數(shù)為：

按照矩形搜索范圍劃分，在修正每個網(wǎng)格節(jié)點時，其搜索面積為一個網(wǎng)格單元的面積，所以平均有N個采樣點用以修正其中心節(jié)點（見圖2），采用該方法最終共可得到三角形面的個數(shù)為2×（columns×rows）。

圖1 標準網(wǎng)格模型劃分

圖2 節(jié)點修正搜索域劃分

2 計算網(wǎng)格模型節(jié)點的修正值

設(shè)在節(jié)點P（r，c）的搜索范圍內(nèi)共搜索到n個采樣點，分別用 Q1， Q2， Q3，…，Qn表示， 它們共同決定了P（r，c）點處的值。設(shè)Q1，Q2，Q3，…，Qn對P（r，c）點的影響系數(shù)分別用ε1，ε2，ε3，…，εn表示， 到P（r，c）點的水平投影距離分別用d1，d2，d3，…，dn表示，

上述n個點共同決定了P（r，c）的值，所以其影響系數(shù)之和為1，即：

設(shè)P（r，c）矩形搜索范圍內(nèi)所有的采樣點與點P（r，c）在水平投影面的距離之和為s，即：

點Qi和P（r，c）在水平面上投影距離用di表示，當di越小，Qi就越接近P（r，c），因而Qi在修正P（r，c）過程中所占的權(quán)重越大，其影響系數(shù)也越大，反之其影響系數(shù)也就越大，即εi∝s－di。

設(shè)：

則式（2）可以表示為：

即：

將式（3）代入式（6）可得：

由式（4）、式（7）可得：

因此，按照εi所確定的權(quán)重對P（r，c）進行高度修正，其修正后的高度值為：

在修正每一個網(wǎng)格節(jié)點時，只需要遍歷數(shù)據(jù)的所有采樣點，找到在待修正節(jié)點的搜索范圍內(nèi)的采樣點即可計算出待修正節(jié)點修正后的值。

3 修正算法優(yōu)化

設(shè)某個數(shù)據(jù)的總采樣點數(shù)為n，經(jīng)過預處理后構(gòu)建了row行、column列的網(wǎng)格模型，平均修正每個節(jié)點需要m個采樣點，則其修正算法時間復雜度為O（row×column×n＋m），通常數(shù)據(jù)采樣點數(shù)在105以上，網(wǎng)格總節(jié)點數(shù)通常也不小于104，如果在單核CPU上以串行方式執(zhí)行此算法，其算法時間復雜度將大于O（109），同時在節(jié)點修正過程中，至少需要進行2m次乘方運算與m次開方運算，因此，其執(zhí)行過程非常耗費計算機計算資源，并且耗時長，因而程序運行效率低。針對以上算法的缺點和不足，筆者提出一種常規(guī)的優(yōu)化方法，具體內(nèi)容如下。

由于要進行建模的數(shù)據(jù)量較大，采樣點較多，所以在所有的采樣點中搜索修正網(wǎng)格所需的點會耗時很多?？紤]到復雜的存儲結(jié)構(gòu)不利于后續(xù)并行優(yōu)化的進行，程序中數(shù)據(jù)在內(nèi)存中采用線性表結(jié)構(gòu)存儲。當數(shù)據(jù)量較大時，線性表的排序算法的時間復雜度也相當大。通過分析地質(zhì)體數(shù)據(jù)的存儲得知，采樣點按照行序依次存儲，據(jù)此為每一個待修正節(jié)點建立一個搜索范圍索引，這樣就能夠大大地減小搜索范圍，進而提高程序運行效率。

考慮到模型的相似度等因素，取Δx＝8x，Δy＝3y，即平均每個節(jié)點由24個采樣點計算所得。針對不規(guī)則數(shù)據(jù)每行的采樣點數(shù)不同的情況，設(shè)原始數(shù)據(jù)平均row行，平均每行采樣點數(shù)為col，程序中動態(tài)建立索引表，每個節(jié)點搜索范圍起始點索引差約為4col。

4 優(yōu)化過程并行架構(gòu)

2007年，NVIDA公司推出CUDA（Compute Unified Device Architecture，統(tǒng)一計算設(shè)備構(gòu)架），在CUDA編程模型［7］中引入了主機與設(shè)備的概念。其中CPU為主機，GPU作為設(shè)備或協(xié)處理器。在主機端主要完成簡單數(shù)據(jù)處理和串行計算，而高度復雜的并行處理計算交由設(shè)備完成。CPU與GPU［8］之間的通訊、數(shù)據(jù)傳送等操作則通過CUDA的API存儲器管理函數(shù)來實現(xiàn)，能夠開辟、釋放和初始化設(shè)備中的存儲空間、實現(xiàn)主機與設(shè)備間的數(shù)據(jù)的傳輸。

圖3 GPU的架構(gòu)與并行數(shù)據(jù)分塊示意圖

當確定了程序中算法的并行部分后，將數(shù)據(jù)拷貝到設(shè)備端。在GPU上運行的函數(shù)稱為kernel（內(nèi)核函數(shù)），其并不是一個完整的程序，而是CUDA程序中的并行執(zhí)行部分。圖3所示為在CUDA程序中主機端和設(shè)備端函數(shù)執(zhí)行的順序。從圖3可以看出，CUDA中的線程被組織成一個二級的層次結(jié)構(gòu)：每個核函數(shù)（kernel）對應(yīng)著一個Grid（柵格），一個Grid（柵格）中包含多個Block（線程塊），而每個Block（線程塊）中則有很多的thread（線程）。一個柵格中不同的Block通過不同的線程塊索引（BlockIdx）來標識，同一個Block中的不同thread通過不同的線程索引（threadIdx）標識，這些索引只能通過核函數(shù)中由設(shè)備專用寄存器提供的只讀型內(nèi)建變量訪問。這樣，在核函數(shù)執(zhí)行時每一個線程都有唯一的一個BlockIdx和一個threadIdx組合與之對應(yīng)，在核函數(shù)中，通過這2個索引確定本線程的具體任務(wù)。

在建模過程中，首先在主機端進行預處理，找出網(wǎng)格模型的邊界值，建立修正網(wǎng)格時所需的索引表并劃分網(wǎng)格。將一個網(wǎng)格模型的修正任務(wù)分配到存在著若干個線程塊的柵格中，線程塊的具體規(guī)模由網(wǎng)格模型大小決定。每一個網(wǎng)格節(jié)點對應(yīng)著線程塊中的一個線程，其能夠安全地訪問采樣數(shù)據(jù)并獨立進行本線程多對應(yīng)網(wǎng)格節(jié)點的修正。由于網(wǎng)格修正算法存在著豐富的數(shù)據(jù)并行性和任務(wù)并行性，因而在GPU上執(zhí)行修正算法能夠獲得很大的性能提升。算法具體流程（見圖4）如下：①讀取原始數(shù)據(jù)；②分析原始數(shù)據(jù)，創(chuàng)建規(guī)則網(wǎng)格模型和索引表；③將原始數(shù)據(jù)和網(wǎng)格數(shù)據(jù)拷貝到設(shè)備存儲器中；④啟動核函數(shù)進行GPU加速，修正網(wǎng)格節(jié)點，同時中間變量存到設(shè)備內(nèi)存的存儲器中；⑤將結(jié)果從設(shè)備內(nèi)存中拷貝到主機內(nèi)存中；⑥劃分三角形網(wǎng)格，剔除不合法的三角形；⑦通過OpenGL繪制網(wǎng)格模型。

圖4 算法主要流程圖

5 試驗結(jié)果與分析

在試驗中，通過處理某工區(qū)地質(zhì)體數(shù)據(jù)并對其進行建模以驗證改進算法。該地質(zhì)體數(shù)據(jù)為不規(guī)則數(shù)據(jù)，在程序中，對數(shù)據(jù)進行預處理后，建立了100×100的規(guī)則矩形網(wǎng)格模型，該網(wǎng)格模型單位矩形規(guī)格為95.157×116.994，直接三角形面片數(shù)為20000個。幾何相似度采用對稱Hausdorff距離進行評估，對稱Hausdorff距離實際上是2個網(wǎng)格模型之間的最大誤差，對稱Hausdorff距離越小，2個網(wǎng)格模型之間的幾何相似度越大。試驗步驟分為3部分：①采用未優(yōu)化的串行算法（簡稱CPU算法）進行處理、計算并繪制網(wǎng)格模型，記錄程序建模所需時間；②采用優(yōu)化后的CPU算法進行處理、計算并繪制網(wǎng)格模型，記錄程序建模所需時間；③采用并行算法（簡稱GPU算法）進行處理、計算并繪制網(wǎng)格模型，記錄程序建模所需時間。CPU串行算法和GPU算法程序運行數(shù)據(jù)對比表如表1所示。由表1可知，與一般串行建模算法相比，由于充分發(fā)掘建模算法的數(shù)據(jù)并行性和算法并行性，利用GPU強大的并行計算能力，使建模速度提高了50～60倍。

表1 CPU串行算法和GPU并行算法程序運行數(shù)據(jù)對比表

在建模精確度方面，雖然利用改進算法構(gòu)建的網(wǎng)格模型與實際地形相比存在一定失真（見圖5），但其失真度能夠控制在很小范圍內(nèi)，因而能夠滿足常規(guī)地質(zhì)建模要求。

圖5 程序運行截圖

6 結(jié)語

提出了一種針對不規(guī)則分布點的快速建模方法，該方法能夠大幅提高建模速度，適用于復雜地質(zhì)體的快速建模。該算法能夠滿足地質(zhì)勘探中對大數(shù)據(jù)量網(wǎng)格快速建模的要求，充分體現(xiàn)了該算法在復雜地質(zhì)體建模中的有效性。

［1］潘懋，方裕，屈洪剛.三維地質(zhì)建模若干基本問題討論 ［J］.地理與地理信息科學，2007，23（3）：75～78.

［2］余杰，呂品，鄭昌文.Delaunay三角網(wǎng)構(gòu)建方法比較研究 ［J］.中國圖像圖形學報，2010，15（8）：28～32.

［3］沅寧君，陳俊，于文茂，等.基于三維地震數(shù)據(jù)的并行希爾伯特黃算法研究 ［J］.計算機工程與設(shè)計，2012，32（6）：89～92.

［4］Jing Ming，Mao Pan，Honggang Qu，et al.GSIS：A 3Dgeological multi－body modeling system from netty cross－sections with topology ［J］.Computers ＆Geosciences，2010，36（6）：756～767.

［5］Hao Zhao，Runcai Bai，Guangwei Liu.3DModeling of Open Pit Based on AutoCAD and Application ［J］.Procedia Earth and Planetary Science，2011，3（1）：258～265.

［6］Wycisk P，Hubert T，Gossel W，et al.High－resolution 3Dspatial modelling of complex geological structures for an environmental risk assessment of abundant mining and industrial megasites ［J］.Computers ＆ Geosciences，2009，35（1）：165～182.

［7］David B K.大規(guī)模并行處理器編程實戰(zhàn) ［M］.趙開勇，汪朝輝，程亦超譯.北京：清華大學出版社，2013.

［8］張舒，褚艷利.GPU高性能運算之CUDA ［M］.北京：中國水利水電出版社，2009.