基于多線程多GPU并行加速的最小二乘逆時偏移算法

柯 璇,石 穎,張 偉,張 振,何 偉

(1.東北石油大學(xué)地球科學(xué)學(xué)院,黑龍江大慶 163318;2.中國石油塔里木油田分公司勘探開發(fā)研究院,新疆庫爾勒,841000;3.中石化石油工程地球物理公司南方分公司,四川成都,610041)

地震偏移成像旨在對地下構(gòu)造的反射信號重新歸位,并根據(jù)地震資料刻畫地下構(gòu)造[1-6]。隨著勘探要求的提升,地震成像的目標(biāo)從地層構(gòu)造描述向地層屬性描述轉(zhuǎn)變[7-10]。將近年來快速發(fā)展的最小二乘逆時偏移算法與反演思想結(jié)合,可用于精確描述地層屬性。YAO等[11]提出了基于矩陣描述的最小二乘逆時偏移算法,并指出該方法在消除低頻噪聲的同時更好地聚焦成像能量。郭振波等[12]提出了最小平方逆時偏移真振幅成像方法,并驗證了該方法在真振幅成像方面的明顯優(yōu)勢。

目前,針對最小二乘逆時偏移算法的研究主要集中于改善算法成像效果、提升算法適用性和計算效率等方面。ZHANG等[13]指出,由于地下介質(zhì)是一種變密度的粘彈性介質(zhì),采用常規(guī)的聲波方程模擬地震波場會導(dǎo)致振幅與實際情況匹配不佳,因此提出了一種新的目標(biāo)函數(shù),降低了振幅的不匹配對最優(yōu)化算法的影響,提升了算法的穩(wěn)定性和數(shù)據(jù)的適應(yīng)性。ZHANG等[14]提出了一種不依賴子波的最小二乘逆時偏移策略,可有效降低由子波不匹配引起的噪聲干擾。李慶洋等[15]提出了去均值歸一化互相關(guān)最小二乘逆時偏移算法,該方法修改了目標(biāo)泛函,利用去均值歸一化算法,降低了算法對子波能量的要求,提升了算法的穩(wěn)定性和可靠性。劉學(xué)建等[16]提出了表面多次波最小二乘逆時偏移算法,并將多次波作為有效信號應(yīng)用于成像算法,增加了成像范圍,雖然初次迭代時產(chǎn)生了串?dāng)_,但隨著迭代次數(shù)的增加,串?dāng)_得以消除。

最小二乘逆時偏移算法的迭代流程計算量大,計算效率的提升對推動最小二乘逆時偏移算法發(fā)展至關(guān)重要。多炮數(shù)據(jù)同時計算是提升最小二乘逆時偏移算法效率的有效途徑之一。BERKHOUT[17]定義了“面炮”偏移概念,先對炮集數(shù)據(jù)合成疊加,再進(jìn)行偏移計算,可有效降低偏移計算量,該思路目前廣泛應(yīng)用于最小二乘逆時偏移算法;平面波靜態(tài)編碼[18]、自適應(yīng)奇異譜分析[19]和頻率選擇編碼[20]等多種方法提升了最小二乘逆時偏移計算效率,也有效壓制了由炮集合成計算產(chǎn)生的串?dāng)_;李闖等[21]推導(dǎo)了平面波最小二乘逆時偏移算法,大幅提升了算法的執(zhí)行效率,改善了最小逆時偏移成像質(zhì)量,分析了多種編碼策略,總結(jié)了多種編碼策略的優(yōu)勢。加快迭代誤差下降速度的方法,可降低計算量,在迭代終止條件不變時,可減少迭代計算次數(shù),以提升計算效率。LIU等[22]針對ZHANG等[13]提出的求取迭代步長參數(shù)的問題,給出了解析步長(analytical step length,ASL)公式,提升了迭代算法的誤差下降效率;預(yù)條件和規(guī)則化方法的應(yīng)用,有效加速了迭代算法的收斂,并帶來了計算效率的提升,提高了深部成像分辨率和保幅性,對于不規(guī)則的地震數(shù)據(jù),該方法有更好的適應(yīng)性[23,24]。GPU加速技術(shù)的發(fā)展,從硬件方面提升了地震數(shù)據(jù)處理方法的計算效率:李博[25]、劉紅偉[26]和SHI等[27]實現(xiàn)了基于GPU加速的逆時偏移算法;石穎等[28,29]將GPU加速技術(shù)應(yīng)用于多次波的預(yù)測和衰減算法中;郭雪豹等[30]采用GPU加速技術(shù)實現(xiàn)了基于頻率衰減的全波形反演方法,這為實現(xiàn)最小二乘逆時偏移算法的高性能計算帶來新的契機(jī)。隨著計算機(jī)技術(shù)的發(fā)展,計算設(shè)備也逐漸升級,如集群設(shè)備中,單個高性能計算節(jié)點通常具備多核中央處理器(central processing unit,CPU)和多GPU,但由于最小二乘逆時偏移算法相對復(fù)雜,需頻繁更新迭代參數(shù),從而導(dǎo)致了最小二乘逆時偏移算法仍缺乏一個相對完整的多GPU加速解決方案。

本文提出了一種多線程多GPU并行加速的最小二乘逆時偏移算法,在GPU加速的最小二乘逆時偏移算法的基礎(chǔ)上,利用CPU的多核架構(gòu),創(chuàng)建多線程協(xié)同操作,調(diào)度多GPU進(jìn)行并行加速計算和迭代參數(shù)的更新,降低數(shù)據(jù)傳輸延遲,大幅提升了計算效率。本文對炮集數(shù)據(jù)分塊切割,在GPU端實行粗粒度并行計算,速度提升接近線性。Marmousi2截斷模型和Marmousi模型的測試驗證了該方法的有效性。

1 時空域最小二乘逆時偏移原理

常密度聲波方程表示如下:

(1)

式中:v(x)為速度場;p(x,t,xs)為波場;f(xs,t)為震源函數(shù);t為時間;xs為震源位置。

假設(shè)速度場v(x)由背景速度場vb(x)和擾動速度場vs(x)疊加而成,即v(x)=vb(x)+vs(x)。由波場疊加原理可知對應(yīng)的波場也由背景波場pb(x,t,xs)和擾動波場ps(x,t,xs)疊加而成,即p(x,t,xs)=pb(x,t,xs)+ps(x,t,xs)。令m(x)=[2vs(x)]/[vb(x)],由波恩近似可得:

(2)

式中:dcal(xg,t,xs)為模擬數(shù)據(jù),xg為檢波點位置。公式(2)的具體推導(dǎo)過程見附錄A。 為了簡化表達(dá),公式(2)可采用矩陣向量形式表示為d=Lm,L為波恩正演算子,d為所有炮的模擬數(shù)據(jù)dcal(xg,t,xs)所構(gòu)成的數(shù)據(jù)向量,m是m(x)的向量表達(dá)形式。

傳統(tǒng)的偏移方法認(rèn)為LT是波恩正演算子的伴隨算子,其表達(dá)式如下:

(3)

式中:pr(x,t,xs)為根據(jù)模擬數(shù)據(jù)所得的檢波點波場,將所有炮的成像結(jié)果進(jìn)行疊加計算,即可獲得最終的成像結(jié)果。為簡化表達(dá),公式(3)也可表示為矩陣向量的形式:m=LTd。

根據(jù)最小二乘逆時偏移算法獲得的擾動模型m來建立最小化目標(biāo)函數(shù):

式中:dobs為觀測數(shù)據(jù)的矢量表示形式。

本文采用CLAERBOUT[31]提出的共軛梯度法,由迭代算法獲得擾動模型m的更新,具體迭代流程如下:

式中:r為數(shù)據(jù)殘差;m為擾動模型,也是迭代計算需要求取的目標(biāo)解;Δm為模型域梯度;Δr為數(shù)據(jù)域共軛梯度;sk,Sk分別為第k次迭代中模型域和數(shù)據(jù)域的更新量;α,β分別為修正sk,Sk的步長參數(shù),〈·〉代表向量的點積運算。為提高計算效率,在迭代過程中,本文參照以下公式對梯度進(jìn)行歸一化補償照明[32]:

(13)

式中:γ為穩(wěn)定性系數(shù);s為當(dāng)前炮數(shù);S為總炮數(shù);tmax為時間方向最大采樣點數(shù)。

2 多GPU加速優(yōu)化方法

GPU加速技術(shù)能夠有效提高并行算法的計算效率[33],已在地震數(shù)據(jù)處理領(lǐng)域取得了較為廣泛的應(yīng)用[25-30]。關(guān)于GPU加速技術(shù)的基本流程不再贅述。本文采用基于CUDA平臺的多GPU加速技術(shù),將GPU存儲器優(yōu)化和多線程多GPU并行加速方法應(yīng)用于最小二乘逆時偏移算法,利用GPU內(nèi)部共享存儲器和寄存器等高速存儲器,降低數(shù)據(jù)訪問延遲,提高計算效率,結(jié)合CPU的多核架構(gòu),分配多CPU線程協(xié)同調(diào)用多GPU進(jìn)行加速計算。

2.1 存儲器優(yōu)化

GPU包含多種存儲器,相較于全局存儲器,共享存儲器和寄存器具有更高的數(shù)據(jù)傳輸帶寬,即更高的讀寫效率。因此,采用共享存儲器和寄存器作為數(shù)據(jù)存儲器協(xié)助計算,可獲得更高的計算效率。

GPU加速計算時,將數(shù)據(jù)的網(wǎng)格點劃分為若干個Block(線程塊),各個Block的GPU線程可執(zhí)行一對一的網(wǎng)格點數(shù)值計算。共享存儲器是各個Block的內(nèi)部存儲器,僅限于同一Block內(nèi)的GPU線程訪問,寄存器則為每個GPU線程的私有存儲器,因此,需合理分配存儲器,才能實現(xiàn)數(shù)據(jù)訪問時的提速。

以(1)式的離散表達(dá)式為例(不考慮震源項):

(14)

圖1 Block劃分及共享存儲器分配示意

2.2 多線程多GPU

目前主流的計算設(shè)備中CPU端具備多核架構(gòu),可同時觸發(fā)多線程作業(yè)。同一計算設(shè)備中配備多個GPU設(shè)備即可支持多GPU并行運算。本文提出了多線程多GPU并行加速最小二乘逆時偏移算法,根據(jù)CPU端多線程機(jī)制,使每個CPU線程負(fù)責(zé)一個GPU的作業(yè)管理和數(shù)據(jù)傳輸,將計算任務(wù)和數(shù)據(jù)分塊,傳輸至GPU端,并以作業(yè)發(fā)送的方式,實現(xiàn)多GPU并行運算,具體情況如圖2所示。

最小二乘逆時偏移算法需進(jìn)行迭代計算,計算量隨迭代次數(shù)線性增加,本文采取的多GPU的并行策略為迭代計算時,先對炮集數(shù)據(jù)進(jìn)行分塊,再分派給各個GPU,彼此獨立地進(jìn)行波場模擬計算。該策略既避免了多GPU間波場數(shù)據(jù)的實時交換,又降低了頻繁的數(shù)據(jù)傳輸引起的計算等待,使得多GPU并行計算結(jié)果逼近線性加速效果。

共軛梯度法迭代時,根據(jù)(7)式和(8)式,對數(shù)據(jù)域維度(炮數(shù)×對應(yīng)炮的道數(shù)×?xí)r間采樣點數(shù))的多個向量進(jìn)行線性運算可求取參數(shù)α和β。編程實現(xiàn)時,如采用BLAS庫運行向量的線性運算,需在GPU端的顯存和CPU端的內(nèi)存之間頻繁地進(jìn)行數(shù)據(jù)傳輸,并且因等待數(shù)據(jù)傳輸而降低算法執(zhí)行的效率。

圖2 多線程多GPU分配示意

向量點積運算滿足(15)式～(17)式所示的性質(zhì),式中A和B分別為兩個向量,對應(yīng)的向量元素分別為a1,a2,…,an,an+1,an+2,…,a2n和b1,b2,…,bn,bn+1,bn+2,…,b2n,該性質(zhì)有利于多GPU并行計算的實現(xiàn),因此可以先分組計算,再對各組結(jié)果求和獲得最終結(jié)果。本文采用CUDA提供的線性代數(shù)程序庫CUBLAS進(jìn)行向量運算,大部分運算均在GPU端執(zhí)行,減少了CPU端內(nèi)存與GPU端顯存之間的數(shù)據(jù)傳輸頻率。

本文采用Pthread接口進(jìn)行CPU端的線程開發(fā),主線程負(fù)責(zé)數(shù)據(jù)同步和狀態(tài)監(jiān)控,利用主線程調(diào)用函數(shù)“Pthread_create()”并創(chuàng)建多個并行子線程后,各子線程負(fù)責(zé)對應(yīng)數(shù)據(jù)塊的讀取,隨后調(diào)用對應(yīng)的GPU進(jìn)行計算。具體實現(xiàn)流程如圖3所示。

圖3 多線程多GPU最小二乘逆時偏移算法流程

主要步驟如下:

第1步,主線程根據(jù)線程個數(shù)進(jìn)行數(shù)據(jù)統(tǒng)計和分塊;

第2步,啟動多線程,各線程根據(jù)任務(wù)分配情況讀取數(shù)據(jù),并將數(shù)據(jù)傳輸至GPU端;

第3步,各線程調(diào)用GPU,根據(jù)(5)式計算梯度;

第4步,設(shè)置多線程阻塞函數(shù)“Pthread_join()”,待所有線程執(zhí)行完畢后調(diào)用CUBLAS庫,對各GPU端的梯度求和并同步;

第5步,各線程調(diào)用GPU,根據(jù)(6)式計算共軛梯度,調(diào)用CUBLAS庫在各GPU端完成共軛梯度法中所需的向量點積計算;

第6步,設(shè)置多線程阻塞函數(shù)“Pthread join()”,待所有線程執(zhí)行完畢后,對各GPU端計算向量點積結(jié)果,并對各GPU所得點積結(jié)果進(jìn)行求和同步,然后根據(jù)(7)式和(8)式計算參數(shù)α和β;

第7步,各線程調(diào)用GPU,并根據(jù)(9)式～(12)式更新迭代結(jié)果及數(shù)據(jù)殘差;

第8步,判斷是否滿足迭代終止條件,如果滿足迭代終止條件,則結(jié)束多線程,輸出數(shù)據(jù);否則重新迭代。

附錄C為C語言編程實現(xiàn)的CPU端多線程作業(yè)觸發(fā)的偽代碼。

3 模型測試

3.1 最小二乘逆時偏移模型測試

3.1.1 模型測試1

本文利用Marmousi2截斷模型進(jìn)行測試,參數(shù)如下:縱、橫向網(wǎng)格點數(shù)分別為296和600,網(wǎng)格間距為15m,雷克子波主頻為16Hz,時間采樣間隔為1.0ms,采樣點數(shù)為6000,設(shè)計20炮震源地表激發(fā),激發(fā)點均勻分布于水平方向1485～7185m,炮間距300m,每炮由199個檢波器接收,檢波器均勻?qū)ΨQ地分布于激發(fā)點兩側(cè),最大偏移距為1485m,檢波器間距15m。

圖4a為準(zhǔn)確速度模型,即正演模型,圖4b為背景速度模型,即偏移算法采用的模型,也是利用準(zhǔn)確速度模型平滑所得的模型,圖4c為擾動模型,可利用準(zhǔn)確速度模型和背景速度模型計算得到,該擾動模型可視為最小二乘逆時偏移的理論解。圖5a為常規(guī)逆時偏移的成像結(jié)果;圖5b為拉普拉斯去噪后的常規(guī)逆時偏移成像結(jié)果;圖5c為最小二乘逆時偏移(50次迭代后)成像結(jié)果。

圖4 模型參數(shù)a 準(zhǔn)確速度模型; b 背景速度模型; c 擾動模型

圖5 逆時偏移處理后得到的成像結(jié)果a 常規(guī)逆時偏移; b 拉普拉斯去噪后的常規(guī)逆時偏移; c 最小二乘逆時偏移(50次迭代后)

對比圖5a,圖5b和圖5c可以看出,相較于常規(guī)逆時偏移,最小二乘逆時偏移能夠獲得分辨率更高的成像結(jié)果,能量更收斂,照明范圍也明顯更廣闊,且振幅與理論值相近,所得結(jié)果具有較為明確的物理意義。

抽取水平方向3km處不同迭代次數(shù)的最小二乘逆時偏移單道數(shù)據(jù)進(jìn)行對比,結(jié)果如圖6所示,隨著迭代次數(shù)的增加,振幅和相位匹配逐漸逼近理論值。

圖6 水平方向3km處單道不同迭代次數(shù)的最小二乘逆時偏移結(jié)果對比a 迭代1次; b 迭代10次; c 迭代50次

圖7所示為數(shù)據(jù)殘差的下降曲線,可以看出,本文方法可對殘差數(shù)據(jù)進(jìn)行有效更新,數(shù)據(jù)殘差隨迭代次數(shù)增加而降低。

圖7 數(shù)據(jù)殘差下降曲線

3.1.2 模型測試2

為進(jìn)一步驗證本文方法的適用性,我們基于Marmousi模型進(jìn)行測試,參數(shù)如下:縱、橫向網(wǎng)格點數(shù)分別為384和122,網(wǎng)格間距為15m,雷克子波主頻為16Hz,時間采樣間隔為1.0ms,采樣點數(shù)為3001,設(shè)計20炮震源地表激發(fā),激發(fā)點均勻分布于水平方向1680～4080m,炮間距120m,每炮由199個檢波器接收,檢波器均勻?qū)ΨQ地分布于激發(fā)點兩側(cè),最大偏移距為1485m,檢波器間距15m。

圖8a為準(zhǔn)確速度模型,即正演模型,圖8b為背景速度模型,即偏移算法采用的模型,也是準(zhǔn)確速度模型平滑所得的模型,圖8c為擾動模型m(x)。

圖8 模型參數(shù)a 準(zhǔn)確速度模型; b 背景速度模型; c 擾動模型

圖9 不同迭代次數(shù)的最小二乘逆時偏移成像結(jié)果a 迭代1次; b 迭代10次; c 迭代50次

圖9a,圖9b和圖9c分別為最小二乘逆時偏移1次、10次、50次迭代后的結(jié)果,可以看出,隨著迭代次數(shù)的增加,成像效果顯著提升。

如圖10和圖11所示,將水平方向2.88km處(圖8,圖9中白色虛線處)和深度方向0.45km處(圖8,圖9中紅色虛線處)理論數(shù)據(jù)和不同迭代次數(shù)的最小二乘逆時偏移結(jié)果進(jìn)行對比,可以看出,隨著迭代次數(shù)的增加,本文方法所得結(jié)果的振幅和相位匹配逐漸逼近理論值,也證明了本文方法對不同參數(shù)模型具普適性。

圖10 深度方向0.45km處單道迭代次數(shù)分別為1次(a),10次(b)和50次(c)的最小二乘逆時偏移結(jié)果對比

圖11 水平方向2.88km處單道迭代次數(shù)分別為1次(a),10次(b)和50次(c)次的最小二乘逆時偏移結(jié)果對比

3.2 最小二乘逆時偏移的多GPU儲存器優(yōu)化方法測試

本文分別采用常規(guī)GPU加速方法(調(diào)用全局存儲器)和存儲器優(yōu)化方法(調(diào)用共享存儲器和寄存器)對模型測試1中的數(shù)據(jù)進(jìn)行1次迭代的最小二乘逆時偏移測試,耗時情況如圖12所示,常規(guī)GPU加速方法耗時約73.7s,存儲器優(yōu)化方法耗時約61.2s。將本文提出的存儲器優(yōu)化方法應(yīng)用于最小二乘逆時偏移算法的波場模擬后,計算效率約提升17%。

圖12 常規(guī)GPU加速方法和存儲器優(yōu)化方法耗時對比

分別使用1～4個GPU進(jìn)行1次迭代最小二乘逆時偏移計算。表1所示為使用不同個數(shù)GPU時,分配給各GPU的炮集數(shù),以1個GPU計算耗時為參考值,計算各個GPU理論上的耗時情況,當(dāng)GPU個數(shù)為3時,分配給GPU2設(shè)備的炮集數(shù)為7,則理論上,GPU2的耗時應(yīng)該為GPU2計算20個炮集數(shù)據(jù)所用時間的7/20倍,以此類推。模型測試1中,參與計算的GPU個數(shù)與耗時情況如圖13所示,黑色實線為實際耗時情況,可看出計算耗時隨GPU個數(shù)增加而降低,多GPU加速效果較為明顯;紅色點劃線為根據(jù)表1中的分配情況,多GPU并行算法達(dá)到理想化完全線性加速時理論上的耗時情況;藍(lán)色虛線所示為各個GPU計算每個炮集數(shù)據(jù)的平均耗時。從圖13中可看出,本文方法的執(zhí)行效率較高,耗時情況接近線性加速,但隨著GPU個數(shù)的增加,算法的實際執(zhí)行效率有所降低。這是由于多GPU并行計算時,需等待最慢的線程完成計算任務(wù),還需要執(zhí)行梯度、更新步長等參數(shù)的同步計算,因此無法完全實現(xiàn)線性加速,從而導(dǎo)致隨著GPU個數(shù)的增加,平均每個炮集數(shù)據(jù)的計算耗時略有增加。

表1 多GPU任務(wù)分配情況

圖13 多GPU實際、理論和平均計算耗時對比(基于模型測試1)

模型測試2的模型尺寸和數(shù)據(jù)量均小于模型測試1的數(shù)據(jù)規(guī)模,參與計算的GPU個數(shù)與耗時情況如圖14所示,我們發(fā)現(xiàn)不同個數(shù)的GPU參與計算時,實際耗時與理論耗時之間的差距會隨著GPU個數(shù)的增加而增大,模型測試2中二者差距的增幅大于模型測試1中的對應(yīng)參數(shù);另外,模型測試2中,平均耗時隨GPU個數(shù)的增加而提升的幅度大于模型測試1中的情況。由此說明,隨著模型尺寸和計算數(shù)據(jù)量的增加,多線程多GPU并行計算方法的加速效率逐漸逼近線性加速,主要原因在于隨著計算量的提升,并行算法的加速作用在整個計算過程中相對增加了,而線程等待和數(shù)據(jù)同步等降低并行效率的計算相對減少了。

圖14 多GPU實際、理論和平均計算耗時對比(基于模型測試2)

4 結(jié)論與認(rèn)識

將GPU加速技術(shù)應(yīng)用于二維時域最小二乘逆時偏移算法,在CPU端采用多線程模式,調(diào)用GPU,實現(xiàn)了多GPU的并行加速計算,迭代計算時,調(diào)用了CUBLAS庫函數(shù)協(xié)助計算,在GPU端實現(xiàn)了數(shù)據(jù)的更新和同步計算,減少了CPU與GPU間的數(shù)據(jù)傳輸,降低了延遲,明顯提升了多線程多GPU最小二乘逆時偏移算法的計算效率。本文還采用了訪問速度更快的共享存儲器和寄存器,進(jìn)一步提升了GPU算法的執(zhí)行效率。

本文提出的多線程多GPU并行加速最小二乘逆時偏移的思路也可應(yīng)用于三維數(shù)據(jù)中,但由于三維空間數(shù)據(jù)量大,在進(jìn)行多GPU加速計算時,建議采用多GPU并行加速策略,即對三維空間數(shù)據(jù)進(jìn)行分塊切割,然后交由各GPU并行運算。多線程CPU支持?jǐn)?shù)據(jù)的并行傳輸,可繼續(xù)充分發(fā)揮并行加速方法的優(yōu)勢。此外,多線程多GPU并行加速方法可普遍應(yīng)用于地震資料的迭代類成像算法以及全波形反演算法,下一步的研究重點應(yīng)集中于本文方法在實際地震資料處理中的應(yīng)用。