王一達 趙長海 李 超 張建磊 晏海華 張威毅(北京航空航天大學計算機學院，北京 100191； 中國石油東方地球物理公司物探技術研究中心，河北涿州 072751)

1 引言

Kirchhoff疊前深度偏移(Kirchhoff Prestack Depth Migration，KPSDM)技術是地震數(shù)據(jù)處理中最重要的深度域成像方法之一。該方法高效靈活，且易于輸出地下共成像點道集，被廣泛應用于成像處理和速度分析。

近年來，隨著勘探地質目標越來越復雜，高密度、全方位野外地震采集非常普遍，數(shù)百TB級數(shù)據(jù)總量的工區(qū)屢見不鮮，地震勘探正邁入PB級數(shù)據(jù)時代，地震資料處理周期越來越長[1]，并行化勢在必行[2]。然而在已知的偏移方法中，KPSDM的高效并行難度較大，其主要挑戰(zhàn)來自執(zhí)行過程中重復讀取地震數(shù)據(jù)和旅行時場數(shù)據(jù)帶來的巨大的輸入/輸出(I/O)開銷[3]。對于當前寬方位三維觀測系統(tǒng)，旅行時場一般能達到數(shù)TB的規(guī)模，每一個地震道對應的旅行時場不同，需要在計算時臨時查詢旅行時表，訪問旅行時的I/O模式非常不規(guī)則，給旅行時數(shù)據(jù)的讀寫優(yōu)化帶來巨大的挑戰(zhàn)。

在集群架構上，為滿足不斷增長的計算需求，當前地震數(shù)據(jù)處理中心單個集群已經(jīng)達到了256～1000個節(jié)點的規(guī)模。隨著可編程圖形處理器(GPU)技術的不斷發(fā)展，集群中的一部分計算節(jié)點會配有若干GPU協(xié)處理器來提升節(jié)點的計算能力。在這種大規(guī)模異構計算集群上，節(jié)點處理能力的一致性、硬件的可靠性假設都會失效。

GPU是一種高度并行的、多線程的眾核處理器。相較于CPU，GPU能夠提供更強的浮點計算能力。此外，根據(jù)Gao等[4]的研究，相比于只有CPU的服務器，配有GPU的服務器在單位能耗下能夠提供更加強大的計算能力。GPU的這種優(yōu)勢使得其在高性能計算領域的應用越來越廣。根據(jù)TOP500的統(tǒng)計數(shù)據(jù)，大約有15%的超級計算機都使用了GPU，并且提供了占TOP500中35%的計算能力[5]。在石油物探領域，通用計算GPU也得到了廣泛的應用。例如Shi等[6]利用GPU對疊前時間偏移程序進行了加速，性能得到了明顯提升。隨著諸如CUDA[7]等GPU編程語言的發(fā)展，使用GPU進行通用計算也變得越來越容易，使得利用GPU加速KPSDM軟件成為可能。然而受限于GPU的一些特性，例如顯存空間有限、CPU-GPU之間數(shù)據(jù)傳輸速度的限制、多線程的執(zhí)行模式等，傳統(tǒng)的深度偏移的任務拆分方法并不能有效發(fā)揮GPU的計算能力。

鑒于以上原因，本文提出了異構計算環(huán)境下的KPSDM混合域并行算法，將傳統(tǒng)的僅從成像空間進行任務拆分的方法，擴展到成像空間、輸入數(shù)據(jù)兩個域，消除了任務間的依賴關系，任務的粒度也更細。同時使用CUDA編程語言對偏移的核心算法進行移植，并根據(jù)GPU的特性實施了一系列優(yōu)化，使得軟件能夠適用于大規(guī)模異構集群系統(tǒng)和復雜的共享運行環(huán)境，達到了運行時容錯、高度可擴展、動態(tài)負載均衡的效果。

2 KPSDM計算復雜度分析

地震資料處理時首先在沿Inline和Crossline方向用一定的間隔對地震工區(qū)進行網(wǎng)格化，每一個網(wǎng)格稱作一個面元，面元的劃分以及孔徑的說明見參考文獻[8]?，F(xiàn)代地震勘探多采用多次覆蓋技術，在一次采集過程中每一個面元將被不同的炮檢組合進行多次記錄，所以三維勘探記錄到一個4維數(shù)據(jù)，4維坐標分別是平面坐標x、y、時間t和炮檢距h。

串行算法邏輯結構如圖1所示，KPSDM的旅行時場是炮(檢)點和地下空間點的函數(shù)，三維情況下旅行時場計算量很大。由于偏移過程中要反復使用旅行時，通常預先計算旅行時場，偏移時再根據(jù)炮檢點位置檢索旅行時。在旅行時場計算之后，對于每一個輸入地震道，首先進行地震道的預處理，隨后分別從線、CMP、輸出道、輸出樣點四個維度循環(huán)進行偏移計算。

假設兩個方向上的偏移孔徑相等，對于r個面元，每個面元中有k個地震道，每個輸出地震道的樣點數(shù)為s，并將橢圓用圓形近似，則每個輸入地震道的樣點要映射到πr2ks個輸出樣點上，輸入地震道數(shù)為n，則算法的時間復雜度是O(nr2ks)。近些年一個勘探工區(qū)采集數(shù)據(jù)已經(jīng)達到幾億地震道，一個面元中的地震道從幾十到數(shù)千個不等，每一道數(shù)據(jù)包含數(shù)千個樣點，偏移孔徑一般在幾百米到上萬米范圍內，每一對樣點的映射包括旅行時的查詢與插值、振幅加權、反假頻濾波、積分求和等計算。一般工區(qū)的旅行時場可以達到數(shù)TB的規(guī)模，在如此龐大的數(shù)據(jù)集內反復查詢旅行時場，數(shù)據(jù)訪問負載極高。

圖1 串行算法邏輯結構

KPSDM并行算法設計的難點主要有以下兩點。①如何利用內存和顯存提升旅行時的檢索速度。偏移的計算量與面元數(shù)成正比，而旅行時數(shù)據(jù)量與面元數(shù)也成正比，內存和顯存的容量分別限制了CPU和GPU能夠處理的最大面元數(shù)量。如果一個地震道的偏移計算量不夠大，那么偏移計算的時間就會小于查詢旅行時的時間，導致軟件的性能受限于數(shù)據(jù)訪問瓶頸，并且會降低處理器的資源利用率。②如何減少地震數(shù)據(jù)重復讀取。成像空間與地震數(shù)據(jù)之間是多對多的關系，即成像空間中的每一個成像塊，都是多塊地震數(shù)據(jù)共同疊加的結果；而每一個地震數(shù)據(jù)塊，都會參與到多個成像塊的計算中。另一方面，成像空間與旅行時是一對一的關系，即成像空間中的每一個成像塊都只需要相應的一個旅行時塊。如果各個任務做同樣的成像空間，對地震數(shù)據(jù)進行拆分，那么每個任務就需要重復讀取同樣的旅行時；反之，如果對成像空間進行拆分，每個任務就需要重復讀取同樣的地震數(shù)據(jù)。

3 現(xiàn)有KPSDM并行算法分析

已有的并行KPSDM算法并不能很好地適應不斷增長的數(shù)據(jù)規(guī)模以及新的集群架構，主要表現(xiàn)在以下三個方面。

3.1 任務的拆分方法

Chang等[9]介紹了一種按Inline方向進行滾動偏移的方法，該方法是建立在“所有節(jié)點的內存之和大于一條測線上所有面元對應的旅行時”這個假設上，該方法在問題域上不具備擴展性，當測線對應的旅行時規(guī)模較大時，通常要求添加更多節(jié)點或者降低追蹤點網(wǎng)格密度、減少測線旅行時對內存的需求。

王華忠等[10]和Rastogi等[11]提出按照深度方向切分成像空間，與Inline方向的切分沒有本質區(qū)別，這種方法唯一的好處是便于過濾孔徑外的地震數(shù)據(jù)，對于小孔徑的偏移作業(yè)可以減少讀取冗余地震數(shù)據(jù)，但實際應用中，小孔徑作業(yè)并不常見。

Panetta等[12，13]采用邊長可變的矩形塊分解成像空間。這種任務拆分方法的優(yōu)點是，任務間耦合度降低，但是旅行時和地震數(shù)據(jù)需要反復從單一存儲上讀取，重復讀取的數(shù)據(jù)量與成像空間的分塊數(shù)成正比，導致KPSDM并行程序的可擴展性與單個節(jié)點的計算能力成反比?？蓴U展性是指在存儲和網(wǎng)絡帶寬既定的情況下，軟件的運行效率隨處理器數(shù)量的增加而線性地提高。該文獻提供的算法受到存儲帶寬的限制，導致單個節(jié)點處理器核數(shù)越多，程序就越容易在節(jié)點數(shù)量比較少的時候達到I/O瓶頸。即使繼續(xù)增加計算節(jié)點，程序性能的提升效果也不明顯。

3.2 任務的調度、分配方案

在偏移階段，部分實現(xiàn)會選擇采用經(jīng)理/雇員(Manager/Worker)架構來組織進程。將一條測線對應的成像面元在所有計算節(jié)點分開后，由Master節(jié)點從共享盤陣讀取該測線輸入范圍內的地震數(shù)據(jù)，并廣播給所有Worker節(jié)點，Worker節(jié)點接收到數(shù)據(jù)后對成像空間進行偏移。廣播的副作用非常明顯，即每個節(jié)點對應的成像空間不同，導致節(jié)點間的計算步調不一致，廣播的同步副作用導致大量CPU時間被浪費。為了克服這些缺陷，Li等[14]提出采用進程組來緩解全局廣播帶來的副作用，但是由于在一個組內仍然需要進行廣播，而且組間需要同步，在復雜的運行環(huán)境下，無法做到負載均衡。

3.3 I/O負載壓力

為降低旅行時對內存空間的消耗，Cunhua等[15]提出了位編碼壓縮算法，利用旅行時在深度方向上連續(xù)性較好的特點，讓每一個深度方向只存放三個原始的旅行時，其余旅行時采用4位二進制存放2階段導數(shù)的值。試算發(fā)現(xiàn)該算法雖然能夠達到3.5倍的壓縮比，但是解壓時間太長，與偏移的時間相當。Alkhalif[16，17]和Teixeira等[18]提出計算旅行時的時候降低網(wǎng)格密度、偏移的時候進行插值的方法，犧牲了旅行時場精度，對偏移質量有影響。并且由于旅行時場占用空間非常大，即使壓縮幾倍，依然有TB級的規(guī)模。盛秀杰等[19]基于谷歌文件系統(tǒng)的設計理念設計并實現(xiàn)了PetroV分布式存儲系統(tǒng)，Li等[20]利用Hadoop的MapReduce編程模型對旅行時的排序進行優(yōu)化，并且將地震數(shù)據(jù)和旅行時數(shù)據(jù)存放在HDFS(Hadoop分布式存儲系統(tǒng))上，從而降低了集中存儲的數(shù)據(jù)訪問壓力。但是這種方法只是將部分數(shù)據(jù)的I/O轉移到了額外的HDFS上，并沒有從根本上解決數(shù)據(jù)重復讀取的問題。

4 混合域并行算法設計與實現(xiàn)方案

并行算法及其實現(xiàn)的主要目標是在大規(guī)模集群系統(tǒng)上高效運行，計算系統(tǒng)的體系結構深刻影響著并行算法的設計，一個典型的地震數(shù)據(jù)處理中心計算系統(tǒng)具有如下特征：

(1)一般處理中心單個集群規(guī)模已經(jīng)達到256個節(jié)點，大型處理中心的集群系統(tǒng)規(guī)模已經(jīng)達到數(shù)千個節(jié)點；

(2)使用集中存儲存放地震數(shù)據(jù)；

(3)計算節(jié)點是普通的雙路x86服務器，至少配置1塊本地盤存放臨時數(shù)據(jù)，內存大小為64～256GB，萬兆以太網(wǎng)或者Infiniband互聯(lián)；

(4)GPU集群系統(tǒng)已經(jīng)是常規(guī)配備，主要用于偏移計算，目前常見的配置是計算節(jié)點配置兩個CPU、兩塊GPU。

異構集群的硬件體系結構可以參見文獻[8]。KPSDM算法設計面向的是數(shù)千異構節(jié)點組成的大規(guī)模集群，任務的耦合度要低，便于運行時容錯。KPSDM并行算法設計關鍵是要處理好三個并行層次：節(jié)點間的進程級并行、CPU內的線程級并行以及GPU內的數(shù)據(jù)級并行。任務拆分要保證并行層次間的數(shù)據(jù)局部性和層次內的可擴展性。還要處理好兩類協(xié)同計算：協(xié)處理器之間、CPU與協(xié)處理器之間。不同處理器間的協(xié)同計算需要盡量保證計算任務之間的異步性與負載均衡。

4.1 并行算法設計

三維地震KPSDM輸入數(shù)據(jù)可以用4維坐標表示，同樣成像空間可以對應到4維空間I(x，y，t，h)。如果從炮檢距維度拆分任務，則任務之間相互獨立，方便進行并行化，但是這樣的拆分方法對一個計算節(jié)點來說任務粒度過大，有三個原因：①作業(yè)中的炮檢距一般只有幾十個，總的任務數(shù)太少；②一個炮檢距成像空間對應的旅行時場一般有數(shù)TB，由于對外存儲的隨機查詢效率很低，因此一個任務的全部旅行時必須能夠放入內存；③在單個節(jié)點內部，為了保證CPU與協(xié)處理器之間計算任務的異步性以及負載均衡，采用CPU與GPU處理相同面元的策略，并將任務對應的地震數(shù)據(jù)在兩類處理器之間動態(tài)分配。這種策略要求一個任務的旅行時能夠常駐GPU顯存中，否則頻繁地在內存與顯存之間傳輸旅行時數(shù)據(jù)會嚴重影響軟件性能。

基于以上原因，進一步把一個炮檢距的成像空間分解為多個面元塊。面元塊對應的旅行時大小與面元數(shù)量是正比例關系，考慮到一般情況下顯存會小于計算節(jié)點的內存，以顯存中能夠容納的旅行時數(shù)據(jù)量作為劃分任務的依據(jù)，確定單個任務包含的面元數(shù)量。

同時，將地震數(shù)據(jù)進一步拆分，把一個炮檢距的地震數(shù)據(jù)拆分為多個數(shù)據(jù)塊。假設數(shù)據(jù)塊大小為d，單個任務的面元塊對應的旅行時量為m，計算節(jié)點物理內存為M，需要滿足：m+d

從成像空間、輸入數(shù)據(jù)兩個域進行拆分之后，一個偏移計算的任務就是炮檢距h、面元塊C和數(shù)據(jù)塊D組成的三元組：Task(hi，Cj，Dk)，每一個任務都是獨立的，而且總的任務數(shù)足夠大，能夠滿足大規(guī)模并行的要求。KPSDM并行算法的流程如圖2所示。在旅行時場的計算完成之后，從炮檢距、面元塊、數(shù)據(jù)塊三個維度進行循環(huán)，每一次循環(huán)就是一個獨立的偏移任務。每個任務又對該面元塊中的所有的輸出樣點循環(huán)進行偏移計算。

圖2 本文的KPSDM并行算法流程

上述混合域任務拆分策略最核心的優(yōu)勢是便于利用本地內存和磁盤優(yōu)化數(shù)據(jù)訪問。上文已經(jīng)介紹過，各個成像面元塊與輸入數(shù)據(jù)塊之間是多對多的關系，因此所有的面元塊都需要同樣的輸入數(shù)據(jù)。如果每個任務獨立讀數(shù)據(jù)，從集中存儲重復讀取的數(shù)據(jù)量與面元塊數(shù)成正比。目前已知的算法通常采用廣播加速重復數(shù)據(jù)的讀取，但廣播會帶來任務間的同步、難以容錯等問題。在本文的并行算法設計中，每一個任務將使用過的旅行時和數(shù)據(jù)塊緩存在本地內存和磁盤，任務獲取數(shù)據(jù)的優(yōu)先級是：本地內存、本地磁盤、遠程節(jié)點和集中存儲。理想情況下，集中存儲的數(shù)據(jù)只需讀一遍，大大降低了集中存儲的I/O壓力。

為便于分析算法的I/O時間，假設每個炮檢距對應的地震數(shù)據(jù)量相同，則理論I/O時間為

(1)

(2)

式中： DS為地震數(shù)據(jù)量；DT為總的旅行時數(shù)據(jù)量；F為炮檢距數(shù)量；α為節(jié)點間的網(wǎng)絡傳輸速率和本地磁盤帶寬的最大值；β為集中存儲的吞吐率；n為作業(yè)節(jié)點數(shù)。

從式(2)可以看出，計算節(jié)點數(shù)越多，本地存儲的聚合帶寬越高，傳輸旅行時與地震數(shù)據(jù)的傳輸效率越高，從而解決了I/O對并行擴展性的限制。

4.2 并行算法實現(xiàn)

對于旅行時計算的并行化實現(xiàn)，在上述任務拆分策略下，旅行時場的各個面元塊的旅行時計算之間沒有耦合性，按照輪詢調度的方式將任務分配給各個計算節(jié)點。計算節(jié)點只需要在計算開始的時候讀取一次速度場數(shù)據(jù)，速度場數(shù)據(jù)一般在幾百MB到幾GB之間，能夠全部放入內存中，可以被多個計算任務重復使用。各個計算節(jié)點完成任務后將結果存放在集中存儲上。由于旅行時計算的并行化策略比較簡單，不再詳細介紹，主要重點介紹偏移部分的并行算法實現(xiàn)策略。

4.2.1 節(jié)點間的并行任務調度

圖3顯示了節(jié)點間的并行任務調度策略，KPSDM實現(xiàn)采用Manager/Worker并行架構，每個計算節(jié)點只啟動一個Worker進程，Worker進程內啟用多線程共同完成一個任務。整個作業(yè)只有一個主節(jié)點上啟動Manager進程，負責向各個Worker分配任務，單個任務單元為4.1介紹的Task(hi，Cj，Dk)。Worker進程收到任務后分別讀取任務對應的地震數(shù)據(jù)并進行計算，然后將成像結果寫回集中存儲。

算法實現(xiàn)的第一個關鍵技術是負載均衡。實際用于地震數(shù)據(jù)處理的集群系統(tǒng)不僅僅是節(jié)點內的異構，還可能存在系統(tǒng)級的異構，即一個集群只有部分節(jié)點配置GPU，計算節(jié)點CPU型號與核數(shù)也可能不一致。KPSDM采用“動態(tài)異步”任務分發(fā)策略保障負載均衡，Worker進程維護一個任務隊列，該隊列只能存放一個任務，一旦任務被取走執(zhí)行，Wor-ker就向主進程請求新的任務，并放入隊列等待執(zhí)行。計算節(jié)點的計算能力越強，請求任務的頻度就越高，從而實現(xiàn)按計算能力分配任務的效果；任務在隊列中等待的同時，可以預先將任務的地震數(shù)據(jù)與旅行時數(shù)據(jù)讀到本地內存，實現(xiàn)等待任務I/O與當前任務計算異步執(zhí)行的效果。

圖3 節(jié)點間的并行任務調度

算法實現(xiàn)的第二個關鍵技術是用本地存儲“擴展”集中存儲的帶寬。每個任務對應的地震數(shù)據(jù)與旅行時數(shù)據(jù)緩存在本地內存和磁盤，服務于其他需要同樣數(shù)據(jù)的任務，同一份數(shù)據(jù)有可能在多個計算節(jié)點上存在，主節(jié)點記錄根據(jù)每個節(jié)點的負載狀況，選擇最空閑的節(jié)點提供數(shù)據(jù)服務。每個計算節(jié)點既作為數(shù)據(jù)的提供者，又作為消費者，發(fā)揮網(wǎng)絡的雙向帶寬，節(jié)點越多I/O聚合帶寬越高，保障了KPSDM程序的可擴展性。

算法實現(xiàn)的第三個關鍵技術是運行時容錯。對于大規(guī)模集群系統(tǒng)來說，計算節(jié)點發(fā)生硬件或者軟件故障的概率會非常高[21]。基于檢查點的容錯一方面I/O開銷太大，與節(jié)點規(guī)模成正比[22]，此外需要重啟整個應用也帶來很差的用戶體驗。為此，開發(fā)了新的并行與分布式編程模型(Geophysical Pa-rallel Programming，GPP)[23]替代MPI，能夠支持應用程序實現(xiàn)運行時任務遷移。GPP運行時系統(tǒng)周期性檢測作業(yè)進程的健康狀況，一旦某進程發(fā)生故障，采用事件機制向所有進程報告發(fā)生故障的進程號。KPSDM主節(jié)點負責捕獲故障事件，將所有分配給故障Worker且沒有執(zhí)行完的任務，重新分配給其他Worker。

算法實現(xiàn)的第四個關鍵技術是慢節(jié)點的處理。慢節(jié)點也是大規(guī)模集群中常見的問題，借鑒Map-Reduce框架的“任務備份”[24，25]策略處理慢節(jié)點可能拖延作業(yè)運行時間的問題。在所有任務分配完畢后，測算已經(jīng)完成任務的平均執(zhí)行時間，假設平均完成時間為t，如果任務運行時間超過λt，向空閑Worker分配該任務作為備份，其中λ的值一般設置為2，太小會過于靈敏，產(chǎn)生過多的備份任務會造成很高的I/O負載，影響整體性能。

4.2.2 節(jié)點內的并行策略

在單個節(jié)點內部，由于有GPU參與計算，需要在GPU和CPU之間進行任務的分配，以保證最大化地發(fā)揮整個計算節(jié)點的計算性能。GPU的計算需要先將需要的數(shù)據(jù)通過PCIE總線(外圍設備互聯(lián)總線)從內存?zhèn)鬏數(shù)斤@存中。相比于GPU的強大的計算能力，CPU-GPU之間的數(shù)據(jù)傳輸帶寬往往成為軟件的性能瓶頸。因此要盡可能降低顯存內數(shù)據(jù)切換的頻率，減少數(shù)據(jù)傳輸?shù)拈_銷。

本文的CPU與GPU之間的并行策略如圖4所示。節(jié)點收到一個成像任務后，CPU和GPU處理相同的面元塊。計算前需要先將該任務對應的旅行時塊加載到內存中，隨后利用CUDA的數(shù)據(jù)傳輸接口，將這個旅行時塊一次性地通過PCIE總線傳輸?shù)紾PU的顯存中。任務對應的地震輸入數(shù)據(jù)在CPU和GPU之間動態(tài)分發(fā)，最后將CPU和GPU的計算結果疊加輸出。這種分配策略的好處是，由于在計算節(jié)點的內存和本地磁盤上緩存著大量的輸入地震數(shù)據(jù)，因此將旅行時傳入顯存之后，充足的輸入數(shù)據(jù)能夠保證GPU在計算足夠長的時間之后再進行旅行時數(shù)據(jù)的切換。另外，在節(jié)點間的任務分配中，使用了一些優(yōu)化措施，使得相同面元塊的偏移任務盡可能地調度到同一個節(jié)點上，這樣即使更換了新的偏移任務，內存和顯存中的旅行時場仍然可能被重復使用，進一步減少了GPU中旅行時數(shù)據(jù)切換的頻率，處理器利用率更高。

在CPU多核之間的并行策略是，每個CPU核啟動一個線程，采用輪詢調度的方式將面元均分給計算線程，保證每個線程分到等量的計算任務，降低線程之間的同步開銷，所有線程共享輸入數(shù)據(jù)，可以提升緩存命中率。微架構方面另外一個重要的優(yōu)化是向量化，生產(chǎn)用的代碼里面存在比較多的條件分支，是向量化的最大阻礙，必須綜合應用代碼各種優(yōu)化方法移除條件分支，將KPSDM核心代碼進行充分向量化之后，整體性能至少提升一倍。

當節(jié)點配備了多個GPU時，各個GPU之間均分面元，與CPU核間的任務劃分方式相同。對于GPU上的計算并行，使用CUDA編程語言實現(xiàn)偏移的核心算法，并且根據(jù)GPU的特性以及CUDA的規(guī)則從以下幾個方面進行了優(yōu)化。

第一個優(yōu)化策略是任務的CUDA線程之間的任務劃分。對于GPU協(xié)處理器，需要在面元的基礎上將任務粒度進一步細分。CUDA中每32個線程被組織成一個線程組(Warp)，每個面元內的輸出地震道由一個Warp線程配合完成，每個線程計算一部分輸出樣點。采用這種GPU任務劃分方法，一方面可以產(chǎn)生足夠多的任務，提高處理器利用率；另一方面，同一個Warp中的線程訪問的地震數(shù)據(jù)位于連續(xù)的128字節(jié)對齊的顯存塊中，訪問模式符合GPU中全局內存(Global Memory)的合并訪問規(guī)則，可以通過一次訪存操作滿足所有線程的讀取請求，如圖5所示。

圖4 CPU-GPU協(xié)同計算

圖5 全局內存的合并訪問

第二個優(yōu)化策略是利用紋理內存(Texture Memory)加速地震數(shù)據(jù)的訪問。紋理內存中的緩存機制能夠提供比全局內存更快的數(shù)據(jù)讀取速度，但是其存儲容量比較小，并不適合對顯存內大塊的旅行時數(shù)據(jù)進行加速；另外，紋理緩存的數(shù)據(jù)是只讀的，也不適合存放成像結果數(shù)據(jù)。地震輸入數(shù)據(jù)在顯存中所占的空間較小，而且在計算過程中要被所有的計算線程只讀訪問，適合利用紋理內存提高訪問速度。

第三個優(yōu)化策略是減少每個線程使用的寄存器數(shù)量，提高處理器的利用率。每個GPU設備只提供數(shù)量有限的寄存器，因此每個線程要求的存儲單元越多，處理器中能夠駐留的計算線程就越少，處理器的利用率也就越低。實現(xiàn)中，由于CUDA函數(shù)的參數(shù)列表中的參數(shù)占用了大量的寄存器，因此將參數(shù)列表中的一部分參數(shù)，尤其是一些在軟件運行過程中不會被改變的指針變量，都存放在常量存儲器中，大幅減少了寄存器的使用數(shù)量。

第四個優(yōu)化策略是利用CUDA流(Stream)隱藏內存到顯存的數(shù)據(jù)傳輸延遲。在偏移計算過程中，盡管旅行時數(shù)據(jù)可以長時間駐留在顯存中，GPU處理器仍然需要等待輸入地震數(shù)據(jù)依次從內存?zhèn)鬏數(shù)斤@存中之后才能執(zhí)行偏移計算過程。CUDA的Stream原理如圖6所示，單個Stream中的傳輸和計算必須按順序執(zhí)行，但是相鄰的不同Stream中的傳輸和計算過程可以并行執(zhí)行。如果計算量足夠大，GPU的計算時間就能夠“覆蓋”數(shù)據(jù)傳輸?shù)臅r間，隱藏數(shù)據(jù)傳輸?shù)难舆t。

圖6 CUDA的Stream隱藏數(shù)據(jù)傳輸?shù)难舆t

在實際生產(chǎn)中發(fā)現(xiàn)GPU也會發(fā)生兩類常見的故障： ①檢測不到GPU； ②雖然檢測到GPU，但運行程序崩潰，主要是由硬件Bug引起，在較新的GPU內比較常見。因此，GPU的容錯非常必要。由于CPU和GPU間的任務耦合性非常低，可以采用一種簡單有效的容錯策略：調用GPU之前，進行一系列簡單的快速測試，如果測試不通過，只用CPU做成像任務。

5 實驗結果與分析

5.1 可擴展性與數(shù)據(jù)吞吐分析

實驗集群有256個計算節(jié)點，每個節(jié)點配置2路8核Intel E5-2670 CPU，128GB內存，2塊500GB的高速HDD(硬盤驅動器)本地盤做RAID0(磁盤陣列)。節(jié)點間50GB/s的Infiniband網(wǎng)絡互連。共享存儲采用GPFS并行文件系統(tǒng)，實測最高并發(fā)吞吐率為4GB/s，超過2GB/s讀延遲會顯著增加。

實驗選用中國東北某工區(qū)數(shù)據(jù)片段，共134GB，約2292萬道，每道有1501個樣點；共79591個面元；炮檢距共73個；總的旅行時場數(shù)據(jù)量為466GB；每一個面元塊有943個面元，對應的旅行時場數(shù)據(jù)量為20GB；數(shù)據(jù)塊大小為51GB。

實驗作業(yè)共生成4250個任務，表1列舉了作業(yè)運行期間分別從共享存儲和本地存儲讀取的數(shù)據(jù)量，作業(yè)使用32～128節(jié)點時，只緩存地震數(shù)據(jù)，旅行時I/O依然在共享存儲；超過128個節(jié)點，共享存儲成為瓶頸，地震數(shù)據(jù)與旅行時場都啟用緩存機制。

表1 運行時數(shù)據(jù)來源分布

從表1可以看出，作業(yè)運行期間共產(chǎn)生了32TB左右的數(shù)據(jù)輸入量，相比原始輸入數(shù)據(jù)，擴大了53倍，未啟用旅行時緩存時，本地存儲分擔了90%的I/O壓力，若將地震數(shù)據(jù)與旅行時場的緩存機制完全打開，可以分流97%的I/O。共享存儲與本地存儲的平均讀速率如圖7所示，作業(yè)節(jié)點數(shù)達到128個節(jié)點時，本地存儲總的讀速率為13.3GB/s，平均每個節(jié)點只需要達到106MB/s的讀速率就可以滿足；共享存儲的讀速率為1.4GB/s，如果再增加作業(yè)節(jié)點，旅行時的讀延遲就會增加，共享存儲就會成為性能瓶頸。打開旅行時場的緩存機制后，KPSDM能夠繼續(xù)加速，最多使用237個節(jié)點時，總的平均讀取速率為19.3GB/s，均攤到每個節(jié)點只有83MB/s。

圖7 共享存儲與本地存儲的平均讀速率

由于分布式的本地存儲提供了極高的聚合帶寬，KPSDM程序的可擴展性非常好(圖8)，甚至出現(xiàn)了超線性加速的情況。出現(xiàn)超線性加速比的主要原因是節(jié)點增多后，數(shù)據(jù)的負載更加均衡，讀取數(shù)據(jù)延遲更低。超過192個節(jié)點之后加速放緩的原因是作業(yè)運行時間變短、作業(yè)收尾時間占比增加。

圖8 節(jié)點并行的執(zhí)行時間與加速比

5.2 節(jié)點內的并行測試與分析

計算節(jié)點共16個物理CPU核，32個超線程，線程級的可擴展性測試結果如圖9所示。當線程數(shù)量從1變化到7時，保持著接近線性的加速比，而從8變化到16時，實際加速比曲線與理想加速比曲線的差距越來越大。造成這種情況的主要原因是每個線程在計算時，讀取不同的旅行時場，各個線程之間可以共享的數(shù)據(jù)非常少，造成三級緩存的命中率隨著線程數(shù)量的增多而降低。從16個線程增加到32個線程時，程序依然在加速，根據(jù)圖中的數(shù)據(jù)可以測算，超線程技術對KPSDM的性能提升了近22%。

圖9 線程級并行的可擴展性

為了驗證CPU與GPU的協(xié)同計算效果，分別使用了NVIDIA的Tesla K40和Tesla K80的GPU測試在單個節(jié)點上CPU與GPU在多種組合策略下的性能加速情況，結果如圖10所示，可以看出：當搭配單個K40的GPU時，軟件性能比僅使用CPU時提高了1.1倍； 當搭配單個K80的GPU時，軟件性能比僅使用CPU時提高了2.14倍； 理想狀態(tài)下，當搭配兩個K80時，軟件性能應該能夠提高4倍以上。然而實際測試顯示，使用兩個K80時程序性能只提高了3.69倍，因為此時CPU-GPU之間的傳輸已經(jīng)成為瓶頸，得不到充足的數(shù)據(jù)供應，處理器的利用率降低，因而影響了加速比。

圖10 CPU-GPU協(xié)同計算加速效果

為了提高CPU-GPU之間的數(shù)據(jù)傳輸效率，NVIDIA公司已經(jīng)提出了一種新的高帶寬且節(jié)能的數(shù)據(jù)傳輸方案NVLink[26]，其傳輸速度能夠達到傳統(tǒng)PCIE3.0速度的5～12倍。這項技術已經(jīng)被應用在NVIDIA新的Pascal架構的GPU中。因此未來如果搭配使用NVLink傳輸方案的GPU，將能夠獲得更好的加速效果。

6 結論

本文提出了一種異構計算環(huán)境下的KPSDM混合域并行算法。該算法從炮檢距、成像空間、輸入數(shù)據(jù)三個域拆分任務，保證了任務之間沒有耦合，便于實現(xiàn)運行時容錯與慢節(jié)點處理；另外，根據(jù)內存或者顯存空間來確定任務大小，保證了任務的數(shù)量足夠多，有利于充分利用大規(guī)模集群的計算資源，實現(xiàn)負載均衡。

對于巨量的地震數(shù)據(jù)與旅行時場帶來的數(shù)據(jù)I/O問題，利用本地內存和磁盤緩存數(shù)據(jù)，降低了共享存儲的壓力。本地內存和磁盤的聚合帶寬為KPSDM的可擴展性提供了強有力的保障。實驗結果顯示，對于實際的地震數(shù)據(jù)資料，并行KPSDM軟件能夠在256節(jié)點的集群上獲得接近線性的加速比，表明該算法具有良好的擴展性。

對于擁有GPU的計算節(jié)點，任務能夠在節(jié)點內部以地震道為粒度進一步拆分，實現(xiàn)CPU與GPU間的動態(tài)負載均衡； 同時根據(jù)GPU以及CUDA編程語言的特性，從數(shù)據(jù)傳輸、數(shù)據(jù)訪問以及處理器利用率等幾個方面對GPU上的核心算法進行優(yōu)化，充分發(fā)揮了GPU的計算能力。實驗結果顯示，搭配使用GPU后，單個節(jié)點的程序性能提升達3.69倍，并且在未來新的GPU架構下，仍有上升空間。

參考文獻

[1] Addair T G，Dodge D A，Walter W R et al.Large-scale seismic signal analysis with Hadoop.Computers & Geosciences，2014，66(5)：145-154.

[2] 羅剛，陳繼紅，孫孝萍等.大規(guī)模異構集群地震作業(yè)調度與資源管理系統(tǒng)的設計與實現(xiàn).石油地球物理勘探，2017，52(增刊2)：200-205.

Luo Gang，Chen Jihong，Sun Xiaoping et al.Heterogeneous cluster scheduling and resource management system for a large number of seismic data processing jobs.OGP，2017，52(S2)：200-205.

[3] 榮駿召，蘆俊，李建峰等.矢量Kirchhoff疊前深度偏移.石油地球物理勘探，2017，52(6)：1170-1176.

Rong Junzhao，Lu Jun，Li Jianfeng et al.Vector pre-stack depth migration based on Kirchhoff integral equation.OGP，2017，52(6)：1170-1176.

[4] Gao Y，Iqbal S，Zhang P et al.Performance and power analysis of high-density multi-GPGPU architectures：A preliminary case study.IEEE 17th International Conference on High Performance Computing and Communications (HPCC)，2015，66-71.

[5] Jack D,Erich S,Horst S et al.Top 500 List.http：//www.top500.org/，2017.

[6] Shi X,Li C,Wang S et al.Computing prestack Kirchhoff time migration on general purpose GPU.Computers & Geosciences，2011，37(10)：1702-1710.

[7] NVIDIA.Cuda C Programming Guide.http：//docs.nvidia.com/cuda/cuda-c-programming-guide/index.html，2013.

[8] 趙長海，羅國安，張旭東等.大規(guī)模異構集群上Kirchhoff疊前時間偏移并行算法.石油地球物理勘探，2016，51(5)：1040-1048.

Zhao Changhai，Luo Guoan，Zhang Xudong et al.Kirchhoff prestack time migration on large heterogeneous computing systems.OGP，2016，51(5)：1040-1048.

[9] Chang H，Van Dyke J P，Solano M et al.3-D prestack Kirchhoff depth migration：From prototype to production in a massively parallel processor environment.Geophysics，1998，63(2)：546-556.

[10] 王華忠，劉少勇，孔祥寧等.大規(guī)模三維地震數(shù)據(jù) Kirchhoff 疊前深度偏移及其并行實現(xiàn).石油地球物理勘探，2012，47(3)：404-410.

Wang Huazhong，Liu Shaoyong，Kong Xiangning et al.3D Kirchhoff PSDM for large-scale seismic data and its parallel implementation strategy.OGP，2012，47(3)：404-410.

[11] Rastogi R，Srivastava A，Khonde K et al.An efficient parallel algorithm：Poststack and prestack Kirchhoff 3D depth migration using flexi-depth iterations.Computers & Geosciences，2015，80(7)：1-8.

[12] Panetta J，de Souza Filho P R P，Da Cungua Filho C A et al.Computational characteristics of production seismic migration and its performance on novel processor architectures.19th Symposium on Computer Architecture and High Performance Computing (SBAC-PAD 2007)， IEEE，2007，11-18.

[13] Panetta J，Teixeira T，de Souza Filho P R P et al.Accelerating time and depth seismic migration by CPU and GPU cooperation.International Journal of Parallel Programming，2012，40(3)：290-312.

[14] Li J，Hei D，Yan L.Partitioning algorithm of 3-D prestack parallel Kirchhoff depth migration for imaging spaces.Eighth International Conference on Grid and Cooperative Computing，2009，276-280.

[15] Cunha C A，Pametta J，Romanelli A et al.Compres-sion of traveltime tables for prestack depth migration.SEG Technical Program Expanded Abstracts，1995，14：180-183.

[16] Alkhalifah T.The many benefits of traveltime compression for 3D prestack Kirchhoff migration.68th EAGE Conference & Exhibition Extended Abstracts，2006.

[17] Alkhalifah T.Efficient traveltime compression for 3D prestack Kirchhoff migration.Geophysical Prospecting，2011，59(1)：1-9.

[18] Teixeira D，Yeh A，Gajawada S.Implementation of Kirchhoff prestack depth migration on GPU.SEG Technical Program Expanded Abstracts，2013，32：3683-3686.

[19] 盛秀杰， 金之鈞， 彭成.PetroV分布式數(shù)據(jù)存儲與分析框架設計.石油地球物理勘探，2017，52(4)：875-883.

Sheng Xiujie，Jin Zhijun，Peng Cheng.PetroV distri-buted-data storage and analytics framework design.OGP，2017，52(4)：875-883.

[20] Li C，Wang Y，Yan H et al.High performance Kirchhoff pre-stack depth migration on Hadoop.Procee-dings of the Symposium on High Performance Computing，Society for Computer Simulation International，2015，158-165.

[21] Schroeder B and Gibson G A.A large-scale study of failures in high-performance computing systems.IEEE Transactions on Dependable and Secure Computing，2010，7(4)：337-350.

[22] Cappello F，Geist A，Gropp W et al.Toward exascale resilience：2014 update.Supercomputing Frontiers and innovations，2014，1(1)：5-28.

[23] 趙長海，晏海華，王宏琳等.面向地震數(shù)據(jù)的并行與分布式編程框架.石油地球物理勘探，2010，45(1)：146-155.

Zhao Changhai，Yan Haihua，Wang Honglin et al.Seismic data processing oriented parallel and distributed programming framework.OGP，2010，45(1)：146-155.

[24] Dean J and Ghemawat S.MapReduce：simplified data processing on large clusters.Communication of the ACM，2008，51(1)：107-113.

[25] Qi C，Cheng L and Zhen X.Improving MapReduce performance using smart speculative execution strategy.IEEE Transactions on Computers，2014，63(4)：954-967.

[26] NVIDIA.NVIDIA NVLink.http：//www.nvidia.com/ object/nvlink.html，2016.