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

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

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

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

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

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

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

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

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

(1)

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

假設速度場v(x)由背景速度場vb(x)和擾動速度場vs(x)疊加而成,即v(x)=vb(x)+vs(x)。由波場疊加原理可知對應的波場也由背景波場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)為模擬數據,xg為檢波點位置。公式(2)的具體推導過程見附錄A。 為了簡化表達,公式(2)可采用矩陣向量形式表示為d=Lm,L為波恩正演算子,d為所有炮的模擬數據dcal(xg,t,xs)所構成的數據向量,m是m(x)的向量表達形式。

傳統的偏移方法認為LT是波恩正演算子的伴隨算子,其表達式如下:

(3)

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

根據最小二乘逆時偏移算法獲得的擾動模型m來建立最小化目標函數:

式中:dobs為觀測數據的矢量表示形式。

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

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

(13)

式中:γ為穩定性系數;s為當前炮數;S為總炮數;tmax為時間方向最大采樣點數。

2 多GPU加速優化方法

GPU加速技術能夠有效提高并行算法的計算效率[33],已在地震數據處理領域取得了較為廣泛的應用[25-30]。關于GPU加速技術的基本流程不再贅述。本文采用基于CUDA平臺的多GPU加速技術,將GPU存儲器優化和多線程多GPU并行加速方法應用于最小二乘逆時偏移算法,利用GPU內部共享存儲器和寄存器等高速存儲器,降低數據訪問延遲,提高計算效率,結合CPU的多核架構,分配多CPU線程協同調用多GPU進行加速計算。

2.1 存儲器優化

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

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

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

(14)

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

2.2 多線程多GPU

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

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

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

圖2 多線程多GPU分配示意

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

本文采用Pthread接口進行CPU端的線程開發,主線程負責數據同步和狀態監控,利用主線程調用函數“Pthread_create()”并創建多個并行子線程后,各子線程負責對應數據塊的讀取,隨后調用對應的GPU進行計算。具體實現流程如圖3所示。

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

主要步驟如下:

第1步,主線程根據線程個數進行數據統計和分塊;

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

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

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

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

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

第7步,各線程調用GPU,并根據(9)式～(12)式更新迭代結果及數據殘差;

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

附錄C為C語言編程實現的CPU端多線程作業觸發的偽代碼。

3 模型測試

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

3.1.1 模型測試1

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

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

圖4 模型參數a 準確速度模型; b 背景速度模型; c 擾動模型

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

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

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

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

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

圖7 數據殘差下降曲線

3.1.2 模型測試2

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

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

圖8 模型參數a 準確速度模型; b 背景速度模型; c 擾動模型

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

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

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

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

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

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

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

圖12 常規GPU加速方法和存儲器優化方法耗時對比

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

表1 多GPU任務分配情況

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

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

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

4 結論與認識

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

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