基于CUDA的測井?dāng)?shù)據(jù)并行處理方法初探

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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

表1 GPU存儲器類型

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

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

3.1 核磁共振測井T2譜反演

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

(1)

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

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

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

XP=A

(2)

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

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

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

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

圖3 T2譜反演處理流程

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

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

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

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

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

4 結(jié) 論

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

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

參考文獻(xiàn):

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

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

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

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

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

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