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

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

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

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

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

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

1 時(shí)空域最小二乘逆時(shí)偏移原理

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

(1)

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

假設(shè)速度場v(x)由背景速度場vb(x)和擾動(dòng)速度場vs(x)疊加而成,即v(x)=vb(x)+vs(x)。由波場疊加原理可知對應(yīng)的波場也由背景波場pb(x,t,xs)和擾動(dòng)波場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為檢波點(diǎn)位置。公式(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ù)所得的檢波點(diǎn)波場,將所有炮的成像結(jié)果進(jìn)行疊加計(jì)算,即可獲得最終的成像結(jié)果。為簡化表達(dá),公式(3)也可表示為矩陣向量的形式:m=LTd。

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

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

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

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

(13)

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

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

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

2.1 存儲(chǔ)器優(yōu)化

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

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

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

(14)

圖1 Block劃分及共享存儲(chǔ)器分配示意

2.2 多線程多GPU

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

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

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

圖2 多線程多GPU分配示意

向量點(diǎn)積運(yùn)算滿足(15)式～(17)式所示的性質(zhì),式中A和B分別為兩個(gè)向量,對應(yīng)的向量元素分別為a1,a2,…,an,an+1,an+2,…,a2n和b1,b2,…,bn,bn+1,bn+2,…,b2n,該性質(zhì)有利于多GPU并行計(jì)算的實(shí)現(xiàn),因此可以先分組計(jì)算,再對各組結(jié)果求和獲得最終結(jié)果。本文采用CUDA提供的線性代數(shù)程序庫CUBLAS進(jìn)行向量運(yùn)算,大部分運(yù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)建多個(gè)并行子線程后,各子線程負(fù)責(zé)對應(yīng)數(shù)據(jù)塊的讀取,隨后調(diào)用對應(yīng)的GPU進(jìn)行計(jì)算。具體實(shí)現(xiàn)流程如圖3所示。

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

主要步驟如下:

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

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

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

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

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

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

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

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

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

3 模型測試

3.1 最小二乘逆時(shí)偏移模型測試

3.1.1 模型測試1

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

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

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

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

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

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

圖6 水平方向3km處單道不同迭代次數(shù)的最小二乘逆時(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)一步驗(yàn)證本文方法的適用性,我們基于Marmousi模型進(jìn)行測試,參數(shù)如下:縱、橫向網(wǎng)格點(diǎn)數(shù)分別為384和122,網(wǎng)格間距為15m,雷克子波主頻為16Hz,時(shí)間采樣間隔為1.0ms,采樣點(diǎn)數(shù)為3001,設(shè)計(jì)20炮震源地表激發(fā),激發(fā)點(diǎn)均勻分布于水平方向1680～4080m,炮間距120m,每炮由199個(gè)檢波器接收,檢波器均勻?qū)ΨQ地分布于激發(fā)點(diǎn)兩側(cè),最大偏移距為1485m,檢波器間距15m。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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