趙 虎 武泗海* 尹 成 唐澤凱 賈 鵬

(①西南石油大學地球科學與技術(shù)學院，四川成都 610500； ②西南石油大學天然氣地質(zhì)四川省重點實驗室，四川成都 610500； ③中國石油東方地球物理勘探公司西南物探分公司，四川成都 610213)

1 引言

地震勘探的發(fā)展與高性能計算息息相關(guān)，勘探新技術(shù)的推廣離不開高性能計算，特別是在疊前深度偏移、逆時偏移和全波形反演等方面，龐大的計算量在一定程度上制約了上述技術(shù)的應用。目前，常規(guī)的CPU加速技術(shù)是基于消息傳遞的MPI編程模型和共享存儲的OpenMP編程模型實現(xiàn)的混合并行[1]。然而，這種加速技術(shù)已經(jīng)無法滿足大規(guī)?？碧接嬎闳蝿盏男枨蟆?/p>

通用計算GPU并行的出現(xiàn)為上述問題提供了良好的解決方案。在地震波場模擬、偏移成像等領(lǐng)域，主流的CUDA和OpenCL異構(gòu)并行編程模型得到廣泛應用。劉紅偉等[2]、李博等[3]實現(xiàn)了在GPU平臺的地震疊前逆時偏移高階有限差分算法。為了將GPU并行用于更大規(guī)模的計算問題，龍桂華等[4]在GPU集群上實現(xiàn)了三維交錯網(wǎng)格有限差分地震模擬；劉守偉等[5]研究了在CPU/GPU集群上進行三維逆時偏移的實現(xiàn)方案。目前GPU并行技術(shù)是解決大規(guī)模計算任務的主要手段之一。

常規(guī)的串行程序無法直接利用GPU進行加速，而只能依據(jù)現(xiàn)有的GPU并行模型重新編寫并行程序。這不僅破壞了原有的串行程序，還增加了并行加速的難度和工作量。OpenACC提供了一種更高效、簡潔的編程模型，且具有良好的硬件兼容性，不僅支持不同廠商的GPU，還支持AMD APU和Xeon Phi等協(xié)處理器。其并行化方式是在原有串行程序的基礎上，通過添加編譯指令實現(xiàn)并行而非重寫代碼。這種并行方式非常類似于OpenMP，編譯器根據(jù)編譯指令對代碼進行并行化，區(qū)別在于OpenMP并行區(qū)域在CPU端執(zhí)行，而OpenACC在GPU端執(zhí)行。

本文首先介紹OpenACC的編程模型和疊前逆時偏移(RTM)的基本原理。引入OpenACC并行技術(shù)，設計實現(xiàn)了逆時偏移的多級并行方案，解決了RTM方法在波場模擬/重建的存儲瓶頸，并通過優(yōu)化數(shù)據(jù)通信，進一步提升計算效率，通過單卡條件下的RTM試驗說明了OpenACC的高效性。最后，通過瓊東南深水坡折帶復雜構(gòu)造模型的逆時偏移成像說明了文中多級并行方案的實用性及高效性。

2 OpenACC編程概述

為了適應更多的不同架構(gòu)設計的GPU設備，OpenACC定義了抽象加速器模型。在抽象加速器模型中，主機/主機內(nèi)存與加速器/加速器內(nèi)存是相互獨立的，兩者不能直接訪問。其執(zhí)行模式為主機端執(zhí)行大部分代碼，而將計算密集型區(qū)域卸載至加速器上執(zhí)行。

目前，大部分加速器支持粗粒度—執(zhí)行單元并行、細粒度—多線程和單指令多數(shù)據(jù)(SIMD)的三級并行。與此對應，OpenACC設計了gang、 worker和vector三個并行執(zhí)行層次(圖1)，其中小方框表示vector通道。對于NVidia設備而言，gang對應于流式處理器(SM或SMX)，vector對應于GPU核心，而worker沒有明確的對應硬件。用CUDA C/C++的術(shù)語來說，gang與block、worker和vector的對應關(guān)系會隨著block的維度而變化。

圖1 OpenACC的3個并行執(zhí)行層次

在OpenACC中，通常利用計算構(gòu)件kernels或parallel并行化計算區(qū)域。kernels適用于自動化并行，生成kernel。而parallel則更適用于需要自主決定線程的布局參數(shù),進行更精細化的并行配置。在kernels并行區(qū)域的變量名保持不變，只是其指代的變量存儲位置發(fā)生變化，數(shù)據(jù)拷貝傳輸已全部由編譯器隱式完成。

由于主機與加速器的存儲相互獨立，因此數(shù)據(jù)管理對于異構(gòu)并行的性能非常關(guān)鍵。在OpenACC中，在編譯指令中添加copyin和copyout子句定義變量的屬性，用于并行區(qū)域數(shù)據(jù)的輸入或輸出。對于更大的共享加速器數(shù)據(jù)并行區(qū)域，可通過添加其他copy( )子句的方式來實現(xiàn)。另外,OpenACC還提供了updata子句來更新局部主機/加速器數(shù)據(jù)。

3 逆時偏移算法分析及并行化

3.1 逆時偏移成像方法

逆時偏移是通過雙程波波動方程對地震資料進行時間上反向外推，并結(jié)合成像條件實現(xiàn)偏移，它避免了上下行波的分離，且不受傾角的限制，并能夠?qū)θ我鈴碗s構(gòu)造進行成像，是目前成像精度最高的地震偏移成像方法。該方法的實現(xiàn)流程主要包括：(1)檢波點波場的反向外推；(2)震源波場的重構(gòu)計算；(3)應用互相關(guān)成像條件，獲得單炮成像結(jié)果；(4)多炮疊加，輸出深度域成像結(jié)果。

本文采用交錯網(wǎng)格高階有限差分算法求解聲波一階速度—應力波動方程[5-7],從而實現(xiàn)檢波點波場外推和震源波場的重構(gòu)計算。

(1)

(2)

式中:p、v分別為地震波場的壓力分量和速度分量；V、ρ分別為介質(zhì)的縱波速度和密度； Δx、Δz、Δt分別為縱、橫向空間采樣間隔和時間采樣間隔；C表示交錯網(wǎng)格差分系數(shù)。

逆時偏移是典型的大計算量、大存儲量的地震數(shù)據(jù)處理方法，其計算瓶頸主要就是波場外推/重構(gòu)計算和成像兩部分[8],這兩部分均具有良好的并行基礎，特別適合于GPU多線程的執(zhí)行模式，計算效率可以得到極大的提高。對于存儲問題，通常采用隨機散射邊界存儲策略,以增加適當?shù)挠嬎懔繛榇鷥r，幾乎不需要增加存儲量，但在成像結(jié)果中會引入頻率較低的噪聲，影響最終成像效果[9]。綜合考慮，本文采用有效存儲邊界策略，在震源正傳過程中采用PML吸收邊界對外部波場進行衰減，由于波場衰減過程的不可逆性，需要同時通過存儲邊界波場外側(cè)一定范圍的每個時刻的外部波場，從而在反推過程中利用存儲的外部波場保證內(nèi)部波場計算的準確性；該策略既有效降低了存儲需求，也避免了在成像結(jié)果引入噪聲。

整個RTM計算由CPU/GPU協(xié)同并行計算實現(xiàn)，CPU負責GPU的控制、數(shù)據(jù)準備等，GPU主要進行耗時大的波場外推及成像計算。最后對互相關(guān)成像結(jié)果采用拉普拉斯濾波[10]去除噪聲。疊前逆時偏移的整體實現(xiàn)步驟為：

①加載震源子波，采用PML吸收邊界，利用高階有限差分進行正向的波場外推，利用GPU對差分計算進行加速，最大計算時間為Tmax,同時保存有效邊界內(nèi)的波場值;

②根據(jù)有效邊界存儲記錄重建震源波場，根據(jù)單炮記錄外推檢波點波場，利用GPU加速計算過程，截止時間為T=0;

③同時讀取步驟②中每個時刻的震源波場和檢波點波場，應用成像條件，利用GPU加速成像計算;

④循環(huán)實現(xiàn)多炮偏移，將偏移結(jié)果進行疊加;

⑤利用拉普拉斯濾波壓制成像結(jié)果中的低頻噪聲;

⑥輸出最終成像結(jié)果。

3.2 RTM算法的OpenACC多級并行實現(xiàn)

波場重構(gòu)/外推是逆時偏移算法的計算“熱點“，本節(jié)以波場迭代計算為例，說明OpenACC實現(xiàn)波場延拓計算GPU的加速的優(yōu)化過程。

波場延拓主要包括速度和壓力的迭代更新兩個獨立的過程，且共享一套模型和場值數(shù)據(jù)。首先在并行區(qū)域入口使用copyin加載模型數(shù)據(jù)，使用create創(chuàng)建加速器端場值數(shù)據(jù)。其中copyin用于需要拷入加速器而無需拷貝回主機的數(shù)據(jù)；而create用于創(chuàng)建僅在加速器中使用的數(shù)據(jù)。在data區(qū)域結(jié)束前，需要利用update更新數(shù)據(jù) (如地震記錄)到主機。#pragma acc kernels用于將循環(huán)計算并行化，生成在加速器端執(zhí)行的kernel，并通過loop independent子句告知編譯器循環(huán)間的數(shù)據(jù)存在獨立性，使得循環(huán)得以安全進行并行化。

通常情況下，逆時偏移模型的尺寸受限于單加速器顯存空間的大小。為了解決單卡GPU顯存不足的難題，采用區(qū)域分解技術(shù)對數(shù)據(jù)進行劃分，并分配給不同的加速器進行計算，可以有效地解決規(guī)模較大模型的逆時偏移的問題。

多卡GPU實現(xiàn)OpenACC并行主要有兩種方案:一是基于共享存儲的OpenMP方案，使用OpenMP啟動多線程，將每個線程關(guān)聯(lián)一個加速器，每個加速器承擔部分計算任務，節(jié)點內(nèi)基于OpenMP實現(xiàn)相鄰加速器數(shù)據(jù)的傳遞和同步，實現(xiàn)線程級的逆時偏移GPU加速；二是基于分布式存儲的MPI方案，利用MPI為每個計算節(jié)點分配一個進程，初始化后，調(diào)用不同節(jié)點上的加速器實現(xiàn)進程級的逆時偏移GPU加速。無論是線程級并行，還是進程級并行，由于區(qū)域分解策略，邊界位置的場值計算需要相鄰加速器的場值，因此必須進行數(shù)據(jù)通信。但同時也會增加數(shù)據(jù)傳輸成本，為了最大程度地降低通信對性能的影響，以基于區(qū)域分解技術(shù)的二維模型為例，分別針對上述兩種多卡并行方案進行優(yōu)化。

3.2.1 基于共享存儲的OpenMP方案優(yōu)化

首先明確OpenMP的特點是線程間的主機內(nèi)存共享，可通過更新主機/設備數(shù)據(jù)即可實現(xiàn)通信(如圖2中的步驟1、2所示)。首先利用主機數(shù)據(jù)更新加速器，計算完成后，經(jīng)由線程同步，再將設備數(shù)據(jù)更新至主機。即可完成一次節(jié)點內(nèi)的通信。

圖2 線程間加速器的數(shù)據(jù)通信

根據(jù)加速器的計算和傳輸引擎相互獨立的特性，利用OpenACC的異步隊列執(zhí)行可以將計算與傳輸重疊執(zhí)行，實現(xiàn)各線程中GPU的通信隱藏，其實現(xiàn)流程如圖3所示。首先對加速器計算區(qū)域進行劃分，依次為C、D、A、B區(qū)域。Stream0執(zhí)行串行流程，H2D和D2H分別表示主機至加速器、加速器至主機的數(shù)據(jù)傳輸。Stream1、Stream2為兩個異步執(zhí)行隊列，前者執(zhí)行計算任務，后者執(zhí)行傳輸任務，藍色箭頭表示執(zhí)行方向，紅色箭頭表示依賴關(guān)系。C區(qū)域的計算與B區(qū)域的數(shù)據(jù)傳輸、D區(qū)域的計算與A區(qū)域的數(shù)據(jù)傳輸均得到重疊執(zhí)行，實現(xiàn)了隱藏傳輸，達到提升計算效率的目的。每個OpenMP線程開始時都要為其綁定一個加速器，async子句為當前kernel指定執(zhí)行流的編號。Stream1執(zhí)行C區(qū)域的計算部分,而Stream2執(zhí)行A區(qū)域的傳輸部分，最后通過wait導語對兩個異步執(zhí)行的流進行同步。根據(jù)PG-Profiler性能分析工具對同步/異步執(zhí)行流的對比結(jié)果(圖4)可以看出，同步執(zhí)行當中，Stream21按順序執(zhí)行通信和計算任務,而異步執(zhí)行中Stream21僅負責初始化的數(shù)據(jù)拷貝，而Stream23負責計算任務，Stream25負責線程間的數(shù)據(jù)通信。異步執(zhí)行使得計算和通信在實現(xiàn)了部分重疊，縮短了用時。

圖3 異步隊列—重疊計算與傳輸流程(a)計算區(qū)域劃分； (b)串行執(zhí)行； (c)異步執(zhí)行

圖4 同步隊列(a)/異步隊列(b)的分析結(jié)果(局部)

3.2.2 基于分布式存儲的MPI方案優(yōu)化

節(jié)點間通過MPI進行通信，將通信數(shù)據(jù)進行適當合并，以此降低通信次數(shù)、提高帶寬利用率。其次，借助于MPI非阻塞式通信等方式來隱藏通信時間。如圖5所示，首先對節(jié)點邊界進行波場計算(步驟1中紅色區(qū)域)，然后計算其余位置波場(步驟2的綠色區(qū)域)，并同時進行節(jié)點邊界場值的MPI非阻塞式通信(步驟2綠色箭頭指示通信方向)，從而隱藏通信延遲，最后執(zhí)行等待非阻塞任務執(zhí)行完畢。

隨著GPU-Direct技術(shù)的發(fā)展，多GPU之間已經(jīng)可實現(xiàn)直接的MPI通信,OpenACC提供host_data導語及use_device子句，可直接訪問設備內(nèi)存的指針，以更簡潔、高效的方式實現(xiàn)加速器間的通信。

為了驗證上述并行方案通信優(yōu)化的結(jié)果，以均勻介質(zhì)模型為例，采用x方向區(qū)域分解策略，固定網(wǎng)格數(shù)Nx=6000，隨著z向網(wǎng)格數(shù)Nz的增加，統(tǒng)計波場迭代300次的平均用時，結(jié)果如圖6所示?？梢钥闯觯涸趩螜C雙卡情況下，OpenMP的執(zhí)行效率高于MPI，考慮OpenMP方案基于共享存儲，不需要額外擴充邊界用于通信，尤其是三維高階差分計算，更降低了對顯存的需求。另外，通信優(yōu)化后，重疊通信與計算時間的效率均有所提升。

圖5 節(jié)點間異步通信優(yōu)化方案示意圖

圖6 并行方案的通信優(yōu)化前后性能對比

4 模型試算

4.1 SEG/EAGE模型

為了說明采用有效邊界存儲策略的OpenACC并行方案在震源波場重構(gòu)方面的準確性，以SEG/EAGE鹽丘模型(圖7)為例進行疊前深度偏移試驗。模型空間網(wǎng)格數(shù)為1290×350，x、y方向網(wǎng)格間距均為12m，采用時間2階、空間14階精度的有限差分完成波場模擬與外推，子波主頻為30Hz，時間采樣間隔為1ms,并記錄t=3500ms的震源正演波場快照、震源重構(gòu)波場快照及檢波點反向外推檢波點快照。

圖7 SEG/EAGE鹽丘速度模型

震源波場計算結(jié)果如圖8a所示，波前面正向傳至3500ms，由于受模型中部高速體鹽丘的影響，波場的傳播特征非常復雜；通過有效存儲邊界內(nèi)的記錄重建3500ms的震源波場(圖8b)，圖8d為兩者的差值，其誤差值量級達到10-4以上，這對于成像結(jié)果的影響是可以接受的。圖8c為根據(jù)地震記錄外推重建的波場。對比圖8a～圖8c中抽取的波場值，結(jié)果如圖9所示，由圖可見正傳波場與重構(gòu)波場的波場值完全一致。通過檢波點外推得到的檢波點的上行波與原始波場也基本一致。

圖8 SEG/EAGE鹽丘模型波場快照對比(a)3500ms的震源正傳波場； (b)有效存儲邊界重構(gòu)3500ms的波場； (c)檢波點外推3500ms的波場； (d)正傳波場(圖8a)與重構(gòu)波場(圖8b)的差值

圖9 正傳震源波場、重構(gòu)震源波場與外推波場的場值對比

在此基礎上，采用單GPU并行方式對疊前逆時偏移進行加速優(yōu)化試驗。偏移總炮數(shù)為200，震源深度為12.5m，均勻布設1290個檢波點，道間距為12m，其余參數(shù)與前文一致。對互相關(guān)成像結(jié)果進行濾波處理后得到鹽丘模型的偏移結(jié)果(圖10)。由圖可見，鹽丘構(gòu)造的成像結(jié)果非常清晰，足以反映復雜的地下構(gòu)造。

圖10 單卡并行模式的鹽丘模型偏移結(jié)果

為了說明OpenACC的優(yōu)異加速性能，分別統(tǒng)計了串行、OpenMP和OpenACC的計算用時速度比(圖11)。由于PGI、Intel、GCC等不同編譯器對程序的優(yōu)化程度不同，可能會導致性能表現(xiàn)出現(xiàn)較大差異。為方便對比，均采用PGI17.4編譯器，并將單線程串行計算的加速比定義為1。OpenMP發(fā)揮了多核計算能力，將加速比提升至5.72倍，而OpenACC借助于GPU實現(xiàn)更大規(guī)模的多線程處理，使得加速比可以達到27.94。

圖11 單GPU并行模式加速效率對比結(jié)果

4.2 實例模型

本文以瓊東南深水坡折帶模型為例，實現(xiàn)優(yōu)化的多級并行逆時偏移。如圖12所示，深水坡折帶模型尺寸為35000m×9000m。逆時偏移中采用高階有限差分進行波場模擬與外推，成像條件為震源能量歸一化[11,12]，數(shù)值模擬參數(shù)分別為：空間采樣間隔Δx=Δz=20m,子波主頻為15Hz,采樣間隔為2ms,接收時間為10s，共采用60炮進行偏移成像，炮點深度為20m,均勻布設1750個檢波點，道間距為20m。

根據(jù)前文所述的兩種并行方案，對比不同條件下單炮偏移的計算成本，為了提高計算成本統(tǒng)計結(jié)果的準確性，采用多炮統(tǒng)計的平均值作為該模式的單炮偏移用時(表1)。

圖12 瓊東南深水坡折帶速度模型

表1 多級并行模式執(zhí)行效率對比

最后，利用OpenACC多級并行方案實現(xiàn)對坡折帶模型的疊前逆時偏移，最終成像結(jié)果如圖13a所示。將其與FFD疊前深度偏移(60炮)結(jié)果(圖13b)進行對比。與FFD[13-15]相比，逆時偏移剖面經(jīng)過拉普拉斯濾波后,有效地壓制了相關(guān)噪聲，成像結(jié)果界面連續(xù)性、保幅性較好，深層成像更加清晰，改善了陡傾構(gòu)造、微小構(gòu)造的成像效果。

本文實驗平臺為Ubuntu Linux 16.04，硬件參數(shù)為： CPU采用Intel(R) Core i7-4790K@4.0GHz，內(nèi)存24G，GPU采用GTX-1050Ti×2,顯存大小為4G/GPU。

圖13 逆時偏移(a)及FFD疊前深度偏移(b)成像結(jié)果對比

5 結(jié)論

本文針對疊前逆時偏移中大計算量、大存儲量的問題，提出并實現(xiàn)了基于OpenACC的RTM多級并行加速方案，并采用有效存儲邊界策略，有效地降低存儲需求。針對本文提出的兩種方案分別進行了通信優(yōu)化，有效地提高了疊前逆時偏移成像方法的計算效率。通過上述研究和分析，得出了以下幾點結(jié)論：

(1)OpenACC是一種簡潔、高效的GPU編程模型，通過添加編譯指令的方式，OpenACC將GPU作為計算協(xié)處理器，對數(shù)據(jù)進行大規(guī)模并行處理。模型計算結(jié)果說明，基于OpenACC的并行方案可以極大地提升疊前逆時偏移方法的計算效率和實用性。

(2)通過本文多級并行方案實現(xiàn)多GPU聯(lián)合計算處理，并結(jié)合有效存儲邊界策略，可以有效地解決大規(guī)模計算中節(jié)點內(nèi)顯存不足的問題。

(3)針對多級并行方案的通信優(yōu)化，一方面，OpenMP方案的異步執(zhí)行的效率高于同步執(zhí)行方案，MPI的非阻塞通信效率高于阻塞通信；另一方面，節(jié)點內(nèi)的OpenMP并行的執(zhí)行效率略高于節(jié)點間MPI并行。

(4)基于OpenACC加速坡折帶模型的逆時偏移成像，取得了良好的效果。另外，與單程波深度偏移結(jié)果相比，逆時偏移在復雜構(gòu)造地區(qū)的成像效果更佳，對局部小構(gòu)造的刻畫更為清晰。