基于OpenACC編程模型的逆時偏移多級并行的設計與優化

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

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

1 引言

地震勘探的發展與高性能計算息息相關，勘探新技術的推廣離不開高性能計算，特別是在疊前深度偏移、逆時偏移和全波形反演等方面，龐大的計算量在一定程度上制約了上述技術的應用。目前，常規的CPU加速技術是基于消息傳遞的MPI編程模型和共享存儲的OpenMP編程模型實現的混合并行[1]。然而，這種加速技術已經無法滿足大規模勘探計算任務的需求。

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

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

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

2 OpenACC編程概述

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

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

圖1 OpenACC的3個并行執行層次

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

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

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

3.1 逆時偏移成像方法

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

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

(1)

(2)

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

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

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

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

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

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

④循環實現多炮偏移，將偏移結果進行疊加;

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

⑥輸出最終成像結果。

3.2 RTM算法的OpenACC多級并行實現

波場重構/外推是逆時偏移算法的計算“熱點“，本節以波場迭代計算為例，說明OpenACC實現波場延拓計算GPU的加速的優化過程。

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

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

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

3.2.1 基于共享存儲的OpenMP方案優化

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

圖2 線程間加速器的數據通信

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

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

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

3.2.2 基于分布式存儲的MPI方案優化

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

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

為了驗證上述并行方案通信優化的結果，以均勻介質模型為例，采用x方向區域分解策略，固定網格數Nx=6000，隨著z向網格數Nz的增加，統計波場迭代300次的平均用時，結果如圖6所示。可以看出：在單機雙卡情況下，OpenMP的執行效率高于MPI，考慮OpenMP方案基于共享存儲，不需要額外擴充邊界用于通信，尤其是三維高階差分計算，更降低了對顯存的需求。另外，通信優化后，重疊通信與計算時間的效率均有所提升。

圖5 節點間異步通信優化方案示意圖

圖6 并行方案的通信優化前后性能對比

4 模型試算

4.1 SEG/EAGE模型

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

圖7 SEG/EAGE鹽丘速度模型

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

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

圖9 正傳震源波場、重構震源波場與外推波場的場值對比

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

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

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

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

4.2 實例模型

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

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

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

表1 多級并行模式執行效率對比

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

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

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

5 結論

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

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

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

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

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