張宮, 李寧, 羅超, 武宏亮, 馮慶付

(1.中國石油勘探開發(fā)研究院, 北京 100083; 2.北京大學(xué)地球與空間科學(xué)學(xué)院, 北京 100871)

0 引 言

隨著測(cè)井技術(shù)的發(fā)展,測(cè)井?dāng)?shù)據(jù)量越來越大,處理方法越來越復(fù)雜,這對(duì)用于處理數(shù)據(jù)的計(jì)算機(jī)和處理軟件提出了新的要求。測(cè)井?dāng)?shù)據(jù)處理軟件以單機(jī)為主,處理程序也都為基于CPU的串行程序,在處理信息量較大或算法較為復(fù)雜的數(shù)據(jù)時(shí),處理能力顯得不足。雖然計(jì)算機(jī)的性能在不斷地提高,但還是很難滿足測(cè)井對(duì)數(shù)據(jù)處理速度的要求,尤其是高端測(cè)井?dāng)?shù)據(jù)處理,有時(shí)處理一口井?dāng)?shù)據(jù)的時(shí)間往往達(dá)到幾十分鐘甚至數(shù)個(gè)小時(shí)。

針對(duì)上述問題,研究了在單機(jī)上利用CUDA(Compute Unified Device Architecture)平臺(tái),基于GPU高效處理測(cè)井?dāng)?shù)據(jù)的方法。CUDA技術(shù)是針對(duì)GPU的C語言環(huán)境之一,利用該技術(shù),能夠方便使用GPU處理復(fù)雜的計(jì)算密集型難題。自CUDA平臺(tái)推出以來,基于CUDA平臺(tái)的高性能并行計(jì)算就被廣泛應(yīng)用于地球物理勘探領(lǐng)域[1]。在地震勘探領(lǐng)域,CUDA技術(shù)已經(jīng)被應(yīng)用于解決地震波數(shù)值模擬[2]、偏移成像[3]等問題;在測(cè)井領(lǐng)域,雙側(cè)向測(cè)井有限元正演算法中CUDA技術(shù)的相關(guān)應(yīng)用已有報(bào)導(dǎo)[4],但測(cè)井?dāng)?shù)據(jù)處理中還未見有相關(guān)應(yīng)用。本文分析了測(cè)井?dāng)?shù)據(jù)處理并行計(jì)算的可行性,嘗試了在測(cè)井?dāng)?shù)據(jù)處理中使用CUDA技術(shù)進(jìn)行加速,并以核磁共振測(cè)井T2譜反演處理為例對(duì)加速效果進(jìn)行了分析。

1 測(cè)井?dāng)?shù)據(jù)處理過程的并行性

一般的測(cè)井?dāng)?shù)據(jù)處理程序包括數(shù)據(jù)加載、數(shù)據(jù)計(jì)算、數(shù)據(jù)寫回3個(gè)步驟。其中數(shù)據(jù)加載和寫回由于需要對(duì)硬盤進(jìn)行讀寫操作,只能作為串行程序執(zhí)行,測(cè)井?dāng)?shù)據(jù)的并行處理集中在數(shù)據(jù)計(jì)算的環(huán)節(jié)。并行數(shù)據(jù)處理分為2種模式,一種是任務(wù)并行性;另一種是數(shù)據(jù)的并行性,只有滿足這2個(gè)條件之一,并行處理才會(huì)發(fā)揮優(yōu)勢(shì)。典型的測(cè)井?dāng)?shù)據(jù)采集方法是以井眼為基礎(chǔ),逐深度點(diǎn)的采集儲(chǔ)層地球物理信息,數(shù)據(jù)處理時(shí)每個(gè)深度點(diǎn)的數(shù)據(jù)可以單獨(dú)處理,這個(gè)特點(diǎn)恰好滿足并行計(jì)算數(shù)據(jù)的并行性要求。

圖1 數(shù)據(jù)依賴示意圖

某些特殊的測(cè)井?dāng)?shù)據(jù)處理時(shí),1個(gè)深度點(diǎn)的計(jì)算結(jié)果,需要從臨近的1個(gè)或多個(gè)點(diǎn)進(jìn)行取值而產(chǎn)生數(shù)據(jù)依賴(見圖1),這將破壞數(shù)據(jù)的可并行性。對(duì)于這種情況,仍舊可以采用迭代重合處理的方法進(jìn)行并行化,即把處理任務(wù)按照結(jié)果個(gè)數(shù)進(jìn)行分割,每個(gè)深度點(diǎn)的結(jié)果作為一個(gè)獨(dú)立的運(yùn)算單元,只是取值時(shí)把當(dāng)前深度的值和附近需要的值一并進(jìn)行讀取。該方法雖然在數(shù)據(jù)讀取的時(shí)候稍有性能的浪費(fèi),但對(duì)于計(jì)算密集型的數(shù)據(jù)處理同樣會(huì)帶來處理效率的提高。

2 基于CUDA平臺(tái)的GPU并行數(shù)據(jù)處理

與CPU不同,GPU中大部分晶體管被用來進(jìn)行計(jì)算,只有少量晶體管被用來進(jìn)行指令控制,因此可以在單個(gè)芯片上集成大量的處理單元(見圖2),從而擁有CPU無法相比的計(jì)算性能。以最新的顯卡K80為例,其擁有高達(dá)4 992個(gè)處理單元,單精度峰值處理達(dá)到8.7萬億次。良好的硬件需要有配套的軟件才能發(fā)揮其優(yōu)勢(shì),CUDA平臺(tái)利用較為通用的C語言擴(kuò)展集,為用戶提供了一系列接口,使在GPU上進(jìn)行通用計(jì)算變得簡單可行。

圖2 CPU與GPU芯片設(shè)計(jì)對(duì)比

與CPU上運(yùn)行的程序設(shè)計(jì)最大的不同就是需要考慮存儲(chǔ)器的選擇,GPU能夠使用的存儲(chǔ)器類型很多(見表1),需要根據(jù)實(shí)際情況進(jìn)行選擇,從而使訪存性能達(dá)到最優(yōu)。

GPU擁有成百上千個(gè)計(jì)算單元,但每個(gè)計(jì)算單元所擁有的寄存器和共享存儲(chǔ)空間卻是有限的,設(shè)計(jì)GPU并行程序時(shí),除了考慮任務(wù)和數(shù)據(jù)的并行性外,還要考慮單個(gè)獨(dú)立的計(jì)算單元需要的存儲(chǔ)空間。例如,對(duì)于核磁共振測(cè)井T2譜反演算法,假如輸入的回波串長度為500個(gè)單精度浮點(diǎn)數(shù),布點(diǎn)個(gè)數(shù)為30個(gè),那么1個(gè)計(jì)算單元需要約60 kB的存儲(chǔ)空間。以GTX765M顯卡為例,每個(gè)線程塊寄存器數(shù)量為65 536,共享存儲(chǔ)容量為48 kB,顯然滿足不了1個(gè)計(jì)算單元60 kB的存儲(chǔ)需求,只能退而求其次,選擇相對(duì)較慢的全局存儲(chǔ)或紋理存儲(chǔ)。

表1 GPU存儲(chǔ)器類型

3 測(cè)井?dāng)?shù)據(jù)處理程序的并行化改進(jìn)

在Windows平臺(tái)上改寫了常規(guī)測(cè)井?dāng)?shù)據(jù)處理程序以及核磁共振T2譜反演程序,并對(duì)處理效率進(jìn)行了測(cè)試。

3.1 核磁共振測(cè)井T2譜反演

T2譜反演是核磁共振測(cè)井?dāng)?shù)據(jù)處理的關(guān)鍵步驟,也是利用核磁共振測(cè)井資料進(jìn)行儲(chǔ)層參數(shù)計(jì)算和油氣水識(shí)別的前提條件。核磁共振測(cè)井獲取的每個(gè)回波信號(hào)實(shí)際都是多種弛豫組分的總體效應(yīng)疊加,可以用多指數(shù)函數(shù)模型表示

(1)

式中,A(t)為t時(shí)刻測(cè)到的回波幅度;T2i為第i種弛豫分量的橫向弛豫時(shí)間;Pi為第i種弛豫分量的零時(shí)刻信號(hào)大小。

每一特征弛豫時(shí)間T2i=αi/ρ,都可以用來表征孔隙尺寸的大小;Pi對(duì)應(yīng)某種特征弛豫尺寸孔隙的孔隙度大小。

在反演中,T2i是預(yù)先設(shè)定一系列的值(T2布點(diǎn)),確定各設(shè)定弛豫T2i和特征弛豫組分m后,結(jié)合回波串A(tj),j=1,…,n可以構(gòu)成一個(gè)超定方程組,其矩陣形式

XP=A

(2)

式中,P=(P1,…,Pm)T;A=(A(t1),A(t2),…,A(tn))T;

由上述方程組計(jì)算得出Pi的過程叫做T2譜反演。一系列預(yù)設(shè)的Ti與對(duì)應(yīng)的Pi構(gòu)成了核磁弛豫分布,它們可以用來表征孔隙大小的分布。T2譜反演算法歸根到底是一個(gè)求解超定線性方程最優(yōu)非負(fù)解的問題。通常這個(gè)方程會(huì)非常大,需要利用矩陣分解或迭代方法進(jìn)行求解,方法有奇異值分解法、模平滑法和迭代法,這幾種算法都存在計(jì)算速度慢的問題[6]。為了測(cè)試,在原有的改進(jìn)型截?cái)嗥娈愔捣纸夥?SVD)反演算法基礎(chǔ)上,利用CUDA平臺(tái)將原程序改寫為GPU并行計(jì)算的T2譜反演算法程序。

3.2 并行算法設(shè)計(jì)及優(yōu)化

核磁共振測(cè)井T2譜反演流程如圖3所示,重要的步驟便是進(jìn)行奇異值分解,如果把奇異值分解算法全部放在GPU上進(jìn)行,就不得不考慮存儲(chǔ)空間的問題。在回波串長度為500,T2布點(diǎn)為30個(gè)的情況下,單個(gè)深度點(diǎn)需要的存儲(chǔ)空間遠(yuǎn)遠(yuǎn)大于單個(gè)線程塊最大寄存器和共享存儲(chǔ)容量,只能考慮全局存儲(chǔ)。

圖3 T2譜反演處理流程

另外,GPU和CPU架構(gòu)差異較大,對(duì)于擁有較大的計(jì)算能力且?guī)掃h(yuǎn)大于CPU,但在數(shù)百個(gè)線程同時(shí)執(zhí)行時(shí),數(shù)據(jù)吞吐方面仍舊會(huì)顯得帶寬不足。對(duì)于T2譜反演算法,除了把計(jì)算任務(wù)移植到GPU上外,還要進(jìn)行訪存模式的優(yōu)化。在程序設(shè)計(jì)完畢后對(duì)程序進(jìn)行了訪存模式檢查,盡量采取合并訪存的方式進(jìn)行數(shù)據(jù)訪問。

3.3 并行結(jié)果及效率分析

為了驗(yàn)證程序的正確性和效率,采用某實(shí)測(cè)井?dāng)?shù)據(jù)進(jìn)行測(cè)試,處理井段長度為1 000 m,采樣間隔0.1 m,處理點(diǎn)數(shù)為10 000個(gè)。每個(gè)深度點(diǎn)采集的回波串個(gè)數(shù)為500個(gè),T2反演布點(diǎn)數(shù)目設(shè)置為30,對(duì)數(shù)均勻分布于0.3～3 000 ms之間。測(cè)試使用的計(jì)算機(jī)CPU為酷睿I74 700MQ,擁有8個(gè)邏輯核心,標(biāo)準(zhǔn)頻率為2.4 GHz;對(duì)比用的GPU為英偉達(dá)GTX765M,擁有768個(gè)處理單元,2 GB的全局存儲(chǔ)。原有在CPU上運(yùn)行的程序處理完成耗費(fèi)時(shí)間131 428 ms,進(jìn)行GPU并行化改進(jìn)后處理耗時(shí)為7 248 ms,加速比達(dá)到18倍左右。

圖4 計(jì)算結(jié)果對(duì)比

圖4中第1道為深度道;第2、3道分別為原有在CPU上運(yùn)行的T2譜反演程序計(jì)算的結(jié)果及移植到CUDA平臺(tái)利用GPU計(jì)算的結(jié)果,其結(jié)果幾乎完全一致;第4、5、6、7道分別是用2個(gè)T2譜計(jì)算得到的幾何平均值、總孔隙度、有效孔隙度和可動(dòng)流體孔隙度,可以看出在個(gè)別地方略有差異外,結(jié)果完全一致。分析發(fā)現(xiàn),計(jì)算結(jié)果的差異是因?yàn)镃PU和GPU浮點(diǎn)數(shù)保留位數(shù)不同,但兩者的計(jì)算結(jié)果在可接受誤差范圍內(nèi)。

4 結(jié) 論

(1) 基于GPU的處理程序具有明顯的效率優(yōu)勢(shì),測(cè)試的T2譜反演程序獲得了約18倍的加速比。說明CUDA的GPU并行計(jì)算可以應(yīng)用于測(cè)井?dāng)?shù)據(jù)處理中,在計(jì)算密集型的測(cè)井?dāng)?shù)據(jù)處理方面可以帶來可觀的性能提升。

(2) 在進(jìn)行GPU并行程序設(shè)計(jì)時(shí),除了考慮處理任務(wù)的并行性外,還要根據(jù)實(shí)際情況進(jìn)行存儲(chǔ)器的選擇及訪存的優(yōu)化,才能使GPU的性能充分發(fā)揮。

參考文獻(xiàn):

[1] 張軍華, 臧勝濤, 單聯(lián)瑜, 等. 高性能計(jì)算的發(fā)展現(xiàn)狀及趨勢(shì) [J]. 石油地球物理勘探, 2010, 45(6): 918-925.

[2] 王守進(jìn), 林年添, 丁仁偉, 等. 利用GPU技術(shù)及分塊策略加速地震波場(chǎng)模擬 [J]. 地球物理學(xué)進(jìn)展, 2014, 29(3): 1292-1297.

[3] 張兵, 趙改善, 黃駿, 等. 地震疊前深度偏移在CUDA平臺(tái)上的實(shí)現(xiàn) [J]. 勘探地球物理進(jìn)展, 2008, 31(6): 0427-0433.

[4] 蘇俊, 尚景濤, 鄒長春. 雙側(cè)向測(cè)井有限元正演算法的GPU并行實(shí)現(xiàn) [C]∥中國地球科學(xué)聯(lián)合學(xué)術(shù)年會(huì), 北京: 2014.

[5] 何宗斌, 張宮, 樊鶴. 一種基于并行計(jì)算技術(shù)提高測(cè)井?dāng)?shù)據(jù)處理速度的方法 [J]. 石油天然氣學(xué)報(bào), 2012, 34(7): 76-79.

[6] 張宮. 核磁共振測(cè)井?dāng)?shù)據(jù)處理方法與軟件開發(fā)研究 [D]. 武漢: 長江大學(xué), 2013.