劉 羽,熊壬浩,肖 熠

(桂林理工大學信息科學與工程學院,廣西桂林541004)

大量實際計算經(jīng)驗表明,電磁法正反演計算量大,非常耗時,計算速度難以滿足當前地質勘查的要求。文獻[1]利用單臺主流PC機對101×50網(wǎng)格、20觀測頻率、21測點的二維航空電磁法數(shù)據(jù)進行正演模擬,耗時4.0～4.5h。文獻[2]采用曙光TC5000A單節(jié)點對34×34×33網(wǎng)格、16觀測頻率、25測點的三維MT數(shù)據(jù)進行反演,運行時間長達47.6h。影響計算速度的主要因素有:①顯式的雅可比矩陣計算耗時太多;②當采用有限單元法求解正演響應時,密集的剖分網(wǎng)格會形成巨大的系數(shù)矩陣導致計算量陡增;③有些反演方法采用自適應算法求取正則因子,導致迭代次數(shù)和總體計算量增加。圍繞電磁法計算速度的提升方法,國內(nèi)外地球物理工作者做了很多算法上的改進[3-5],特別是在并行計算方面,進行了很多探索。2006年,譚捍東等[6]基于頻率劃分,采用MPI平臺實現(xiàn)了大地電磁三維正演的并行計算;同年,劉羽等[7-8]基于頻率和拉格朗日乘子劃分,在并行計算機(PVM)平臺上,實現(xiàn)了二維MT Occam并行反演;陳露軍[9]實現(xiàn)了基于MPI的井地三維電磁數(shù)值模擬;2009年,SIRIPUNVARAPORN等[10]利用頻率劃分方式實現(xiàn)了包含垂直磁場轉換函數(shù)的三維數(shù)據(jù)空間并行反演;2010—2012年,胡祥云等將并行計算應用于二維大地電磁正演,然后擴展到三維數(shù)據(jù)空間反演,取得了較好的并行效果[2,11]。

前述并行計算成果均采用了消息傳遞模型和基于頻率劃分的大粒度并行方法,當觀測頻率數(shù)較少時,此種方式不具備可擴展并行性,計算性能難以通過增加節(jié)點數(shù)得到相應提高。例如,文獻[6]在4節(jié)點集群上獲得2.42的加速比(2D prism模型);文獻[2]和文獻[8]在8節(jié)點集群上分別獲得5.69和5.65的加速比;文獻[10]在8節(jié)點和16節(jié)點集群上分別獲得4.5和7.3的加速比(采用PCG求解器)。可見,要想獲得更高計算性能,僅靠擴大集群規(guī)模很難奏效。

本文將傳統(tǒng)的CPU集群擴展為CPU/GPU異構集群,以MT Occam反演為例,基于不同粒度并行分量的挖掘,實現(xiàn)計算任務的不同層次映射,從而在大粒度并行的基礎上進一步提升計算性能。首先給出MT Occam反演的理論方法,然后分析反演方法的并行特征和現(xiàn)有并行方法存在的性能瓶頸,并針對性地提出混合并行方案,利用2種規(guī)模的5個理論模型進行了反演測試,給出了方案的實際效果,最后討論了方案的適用性并給出結論。

1 反演方法

MT Occam反演目標函數(shù)可表示為[12-13]:

(1)

式中:μ為拉格朗日乘子;d為觀測數(shù)據(jù)向量;m為模型參數(shù)向量;F(m)為模型正演響應;W為歸一化矩陣;‖?m‖2為模型粗糙度,對二維構造,同時考慮了模型的橫向和垂向光滑問題,?=[?y,?z]T,其中?y和?z分別表示橫向和縱向相鄰的模型參量的粗糙度矩陣;X2為d與F(m)的擬合差,表示為:

(2)

在迭代尋優(yōu)中,需要在每一步使用局部線性化[2]。設當前模型為mk,則mk+1的響應為:

(3)

式中:Δm=mk+1-mk為模型的修改量;Jk=(?F/?m)k為模型mk處的雅可比矩陣。將(3)式代入(1)式對m求導置零,可得到以μ為變量的模型迭代表達式:

(4)

式中:dk=d-F(mk)+Jkmk。μ值通常使用掃描法或者一維搜索法求得,其取值應使得殘差平方和‖Wd-WF[mk+1(μ)]‖2最小。

MT Occam反演是迭代算法,通過不斷修改給定的初始模型,使其逐步趨近真實模型。每次迭代通過綜合判斷擬合差和粗糙度來找出最優(yōu)解,滿足一定條件時,輸出模型并結束計算。

2 并行性分析

如前所述,反演迭代計算之間嚴格順序依賴關系,無法直接并行實現(xiàn),因此并行計算只能在一個迭代步內(nèi)進行。一個迭代步包含模型正演響應F(m)、雅可比矩陣Jk、叉積矩陣、楚列斯基(Cholesky)分解、拉格朗日乘子μ等主要計算項,不同計算項通常具有不同的并行特征,因而需要采用不同的并行策略。

2.1 粗粒度及細粒度并行性

MT Occam反演的一個迭代步內(nèi),含有粗粒度和細粒度并行性。非常耗時的雅可比矩陣計算和運行次數(shù)極多的模型正演響應計算都是針對頻率分別進行的,各頻率計算之間無數(shù)據(jù)耦合,容易基于頻率劃分實現(xiàn)大粒度并行。當采用掃描方式求取拉格朗日乘子時,μ值點預先按固定間隔選取,其計算相互獨立,也可以實現(xiàn)大粒度分解。μ值一維搜索方式下,采用黃金分割法尋找凹區(qū)間,新μ值點的確定依賴前點,屬于迭代方式,無法直接并行,但可以在正演響應計算部分以分頻方式實現(xiàn)并行計算。

細粒度并行性分布于反演過程中不同的計算項。例如,一個頻率的模型正演響應計算,需要采用高斯消元法求解形如ax=b的大型線性方程組,其中a為有限單元剛度陣,b為常數(shù)項,x為待求向量,消元過程中的行操作可以分解,形成極細粒度計算任務,這些計算可以通過解耦實現(xiàn)并行處理。雅可比矩陣計算也具有同樣的特征。此外,每次迭代中,需要完成多個矩陣叉乘運算和多次Chelosky分解,這些部分也都具有極細粒度的并行性。

2.2 傳統(tǒng)并行方式的性能瓶頸及解決方案

基于頻率和μ值掃描點的分解具有天然的并行性,粒度也可以通過任務組合加以調(diào)整,適合消息傳遞模式下的大粒度并行,是電磁法正反演傳統(tǒng)并行計算的主要分解方式,文獻[2]、文獻[6]和文獻[10]都采用了頻率劃分策略,文獻[8]則結合采用了頻率劃分和μ值掃描點劃分。但一般來說,電磁法勘探中的觀測頻率數(shù)有限,通常不會超過40個,反演中μ乘子掃描點數(shù)一般也不會超過30個,因此,基于頻率和μ值掃描點的大粒度分解方式并行度是受限的[14]。以f表示程序中的串行部分,T1為串行執(zhí)行時間,假定并行工作由N份大小固定的工作tu組成,每份工作可以與其它工作并行執(zhí)行,但不再繼續(xù)細分(一個包含N個觀測頻率的模型正演及偏導數(shù)計算就具有此種特征),則在一個含有p個處理機的并行機上執(zhí)行程序的時間可以表示為:

(5)

公式(5)中N實際為并行度,當P超過這一并行度時,并行運行時間便不再減少。這暴露出MT Occam反演傳統(tǒng)并行方式的性能瓶頸:

1) 當并行度較低時,增加處理機數(shù)沒有意義。

2) 由于并行度并非遠大于處理機數(shù),很難保證所有處理機能幾乎平均地分配計算任務,這容易導致短板效應,降低效率。

由此可見,在消息傳遞模式下,MT Occam反演的并行計算性能由于并行度受限而難以提高。

要繼續(xù)減小并行計算時間Tp,必須細分tu,即充分利用算法中的細粒度并行性。但不同并行特征的計算只有適應計算機集群的體系結構才能獲得最佳的計算效率[15]。在消息傳遞模式下,tu的細分沒有意義,因為計算每個細粒度任務需要進行繁雜的通信,每一次的通信都包含啟動時間和傳輸時間,因此計算/通信比會很小,程序的運行效率非常低下,極端情況下甚至會帶來災難性后果。因此,細粒度并行不適合消息傳遞模式,需要運行在共享內(nèi)存架構上,這需要對傳統(tǒng)結構進行擴展。

由于對細粒度并行的良好支持,GPU的通用計算技術(GPGPU)發(fā)展極為迅速并得到廣泛應用[16-17]。CPU和GPU能夠容易地構成混合異構系統(tǒng),利用傳統(tǒng)消息傳遞模型、OpenMP編譯制導和最新GPU編程接口CUDA可以實現(xiàn)混合并行算法,從而靈活處理具有不同并行特征的計算問題。由于MT Occam反演中存在的層次并行特征,本文采用異構混合并行方案來解決反演中不同粒度的并行計算,利用GPU與CPU的協(xié)同工作,實現(xiàn)系統(tǒng)整體計算能力的最大化。

3 混合并行方案設計

圖1給出了混合并行模型框架,硬件自上向下分別是集群、節(jié)點和GPU計算設備,對應的軟件環(huán)境為MPI,OpenMP和CUDA。

上層分布式并行模型MPI用于實現(xiàn)計算的節(jié)點間并行,本文采用全局靜態(tài)通信,在映射過程中由主節(jié)點進程(主進程)負責對子節(jié)點進程(子進程)進行主動發(fā)送和接收,子進程被動接收和發(fā)送,在確認主進程接收前持續(xù)阻塞狀態(tài)。圖2給出了MPI層主進程的通信過程的實現(xiàn)代碼,其中nproc為子進程數(shù)。

OpenMP層對上層發(fā)送的任務繼續(xù)細分,并對GPU進行管理。反演的主要計算由GPU完成,節(jié)點CPU只負責少量的串行計算,因此沒有設計利用CPU多核的算法來完成GPU之外的工作,加入OpenMP層的主要目的是為將來預留GPU擴展機制。

CUDA層負責完成有限單元法計算、矩陣叉乘、Chelosky分解等具有極細粒度重復計算特征的部分。主機端需要將計算數(shù)據(jù)拷貝至設備端顯存,調(diào)用Kernel函數(shù)通過多線程方式完成并行處理,CPU最后將結果拷貝回內(nèi)存。

圖3為MT Occam CPU/GPU并行反演簡化流程:主進程讀入數(shù)據(jù)、對頻率分組并分發(fā)數(shù)據(jù)到各子進程;子進程接收任務,調(diào)用GPU完成模型正演和雅可比矩陣計算;主進程接收子進程計算結果并對模型正演響應和雅可比矩陣進行匯總,調(diào)用GPU矩陣算法計算叉積矩陣,針對不同的μ值,調(diào)用GPU Cholesky分解算法計算相應模型并發(fā)送給子進程;子進程接收模型數(shù)據(jù),調(diào)用GPU計算正演響應;主進程接收模型正演響應數(shù)據(jù),計算其擬合差和粗糙度,評判模型,確定是進入下次迭代還是輸出最終模型。

圖1 CPU/GPU混合并行模型框架

圖2 MPI層主進程通信代碼

計算流程中正演響應和雅可比矩陣利用有限單元法計算,其線性方程組通過高斯消元法求解,此部分計算占比最大,是提高計算效率的核心部分。本文針對其系數(shù)矩陣的二維等帶寬壓縮存儲結構和MT Occam反演的特點,設計了高斯消元GPU算法。算法以一個線性方程組的消元過程(一個工作三角形)為計算任務映射到CUDA模型的BLOCK線程塊,然后由CUDA二次映射到流處理器上完成,總的映射任務數(shù)與節(jié)點計算頻率數(shù)和GPU個數(shù)相關[18]。

叉積矩陣的求取和Cholesky分解計算量相對較小,不再分發(fā)到計算節(jié)點,而由主進程調(diào)用GPU單獨完成。考慮到NVIDIA公司提供的函數(shù)庫已經(jīng)包含了部分高效算法,本文直接調(diào)用子程序cublasSsyrk_v2計算叉積矩陣,調(diào)用設備端函數(shù)cula_device_SPOSV實現(xiàn)Cholesky分解。

圖3 MT Occam CPU/GPU并行反演簡化流程

4 理論模型試算

本文采用的異構集群配置4個節(jié)點,每個節(jié)點配置1個物理4核i5-3470 CPU,8GB DDR3雙通道內(nèi)存,1個具有2GB顯存的GTX 680 GPU(1536個流處理單元),節(jié)點間利用千兆局域網(wǎng)絡連接。操作系統(tǒng)為CentOS 6.4,開發(fā)工具采用CDK,包含支持MPI,OpenMP和CUDA的C/C++和Fortran編譯器,同時提供了性能分析工具和調(diào)試工具。

4.1 反演實驗模型

為了評估算法的實際運行效率,本文實驗過程中設計了5個模型,在所構建的異構集群上進行試算。每個模型均有2種采樣規(guī)模(即類型1和類型2),限于文章篇幅,這里只展示2個典型模型及其反演結果,即模型1(類型1)、模型1(類型2)、模型2(類型1)和模型2(類型2)。2種采樣規(guī)模的4項主要參數(shù)見表1,表中Nr為測點數(shù),Nm為模型參量,Ne為有限單元網(wǎng)格數(shù),Nf為觀測頻率數(shù)。最大反演次數(shù)定義為20。

表1 理論模型參數(shù)

模型1(圖4)設計為1000Ω·m基底,中部為相同電阻率的地壘,截面尺寸5.0km×12.5km,上部為100Ω·m均勻層。模型2(圖5)設計為100Ω·m均勻半空間有2個不同電阻率棱鏡體,上部為高阻層。

圖4 反演模型1

圖5 反演模型2

類型1模型測點數(shù)為21,模型參量為889,有限單元網(wǎng)格數(shù)為103×35,觀測頻率數(shù)為20。類型2模型加密了各項參數(shù),4項參數(shù)分別增加至41,1549,193×35和36。反演初始模型定義為100Ω·m均勻半空間,理論模型數(shù)據(jù)均加入5%的隨機噪聲。

4.2 算法測試結果

對類型1模型,分別用串行和并行算法計算5個模型的反演。定義偏差值為本次迭代串行算法所計算出的模型擬合差和粗糙度與相對應并行算法計算值之間差值的絕對值。圖6給出了模型2反演過程中擬合差和粗造度的偏差值走勢曲線,偏差值均不大,說明并行計算過程相對穩(wěn)定。

表2給出了類型1和類型2模型串行算法和并行算法的迭代次數(shù)、計算時間和加速比的統(tǒng)計結果。結果表明,規(guī)模較大的模型(類型2)計算時間明顯增多,但加速比也更大。個別模型出現(xiàn)正演次數(shù)的差異,分析原因是CPU與GPU處理器在浮點運算精度上的細微差異所致,最終反演結果沒有受到影響。圖7至圖10展示了模型1和模型2兩種規(guī)模下并行反演的結果(黑色實線為理論模型結構輪廓),對模型1來說,兩種規(guī)模反演結果無太大差別,但從模型2的反演結果看,類型2與理論模型明顯更為接近。

本混合并行方案反演結果正確,與純消息傳遞方案相比,加速比有較大提升,說明方案合理可行。不同模型加速比有一定差異,這取決于模型規(guī)模、復雜度和映射到不同層次計算量的變化?？傮w而言,計算密度大的模型(如類型2)有著更大的加速比。在本文4節(jié)點環(huán)境下,類型1模型平均加速比為6.11,最高加速比為7.11,類型2模型平均加速比為16.46,最高加速比達23.49。

圖6 模型2(類型1)誤差曲線

表2 理論模型反演統(tǒng)計表

圖7 模型1(類型1)反演結果

圖8 模型1(類型2)反演結果

圖9 模型2(類型1)反演結果

圖10 模型2(類型2)反演結果

5 方案適用性分析

影響本文方案性能的因素主要源于兩方面:①通信延遲。本文測試結果表明,混合模型各層均會造成一定的通信延遲。由于OpenMP層和CUDA層是共享內(nèi)存方式,其延遲要遠小于消息傳遞方式,因此,MPI層通信是主要的瓶頸。消息傳遞模型要求計算任務粒度盡可能大以提高通信效率,當模型較小特別是頻點較少時(如類型1模型),本方案上層的任務粒度就相對較小,盡管在下層采用了GPU處理,整體計算性能仍會受到影響。②GPU計算并行度。方案中,反演過程的大量計算最終將映射到GPU的流處理器上完成,細粒度并行度便成為影響性能的因素。GPU計算要求一次將盡可能多的數(shù)據(jù)投入到運算中,當并行度不夠時,流處理器(SP)資源不能得到充分利用,計算效率將會降低。由于上述2方面的影響,本文混合并行方案更適合于規(guī)模較大的模型(多測點、多頻率、大密度網(wǎng)格)反演。

此外,系統(tǒng)存儲性能制約了本文集群的適用范圍。由于MT Occam反演對存儲空間的需求較大,小規(guī)模集群通常適用于二維模型反演,測試結果表明本文設計的不同規(guī)模的多個二維模型反演完全可行。但就三維反演而言,由于有限單元節(jié)點數(shù)遠遠大于二維反演,其需要占用的存儲空間是小型集群存儲資源難以滿足的,因此,如果要用于小規(guī)模三維模型反演,節(jié)點內(nèi)存和GPU顯存都必須增加,并對存儲進行優(yōu)化。

6 結論

MT Occam反演中傳統(tǒng)的并行方法不具擴放性,計算性能的提升難以利用擴大集群規(guī)模的方式來實現(xiàn)。采用CPU/GPU多層次混合并行架構,通過挖掘細粒度并行分量,可以很好地彌補這一不足。線性方程組Gauss消元求解是反演過程中最為核心的計算,在GPU上實現(xiàn)其細粒度并行,能使極為耗時的雅可比矩陣和模型正演響應計算都得到加速,是提升反演速度的關鍵。矩陣叉乘、Cholesky分解也具有良好的細粒度并行性,并可直接利用CUDA并行函數(shù)庫高效計算,這進一步提高了算法的綜合計算性能。

本方案在4節(jié)點的小型集群上對較大模型(類型2)反演獲得了平均16.46的加速比,最高達23.49,相對于純消息傳遞方式有較大提高。較小模型(類型1)的反演受MPI層計算粒度和GPU層并行度的影響,平均加速比為6.11,最高為7.11,效率相對要低,但仍高于純消息傳遞方式?；旌喜⑿心Ｊ较峦ㄐ叛舆t是需要關注的問題,在我們的實驗中,三層模型通信耗時最高達總計算時間的30%左右,說明數(shù)據(jù)通信已經(jīng)成為本方案效率的主要制約因素,需要進行更多的優(yōu)化,特別需要考慮計算與通信的重疊(overlap)。

致謝:感謝MT Occam 反演作者Steven Constable先生提供的串行源代碼。

參 考 文 獻

[1] 李小康.基于MPI的頻率域航空電磁法有限元二維正演并行計算研究[D].北京:中國地質大學,2011

LI X K.A MPI based parallel calculation on two dimensional finite element modeling of AEM[D].Beijing:China University of Geosciences,2011

[2] 胡祥云,李焱,楊文采,等.大地電磁三維數(shù)據(jù)空間反演并行算法研究[J].地球物理學報,2012,55(12):3969-3978

HU X Y,LI Y,YANG W C,et al.Three-dimensional magnetotelluric parallel inversion algorithm using data space method[J].Chinese Journal of Geophysics,2012,55(12):3969-3978

[3] 吳小平,徐果明.大地電磁數(shù)據(jù)的Occam反演改進[J].地球物理學報,1998,41(4):547-554

WU X P,XU G M.Improvement of Occam’s inversion for MT data[J].Chinese Journal of Geophysics,1998,41(4):547-554

[4] SIRIPUNVARAPORN W,EGBERT G.An efficient data-subspace inversion method for 2-D magnetotelluric data[J].Geophysics,2000,65(3):791-803

[5] 陳小斌.大地電磁正反演新算法研究及資料處理與解釋的可視化集成系統(tǒng)開發(fā)[D].北京:中國地震局地質研究所,2003

CHEN X B.New forward and inversion algorithms and a visual integrated system for MT data[D].Beijing:Institute of Geology,China Seismological Bureau,2003

[6] TAN H D,TONG T,LIN C H.The parallel 3D magnetotelluric forward modeling algorithm[J].Applied Geophysics,2006,3(4):197-202

[7] 劉羽,王家映,孟永良.基于PC機群的大地電磁Occam反演并行計算研究[J].石油物探,2006,45(3):311-315

LIU Y,WANG J Y,MENG Y L.PC cluster based magnetotelluric 2-D Occam’s inversion parallel calculation[J].Geophysical Prospecting for Petroleum,2006,45(3):311-315

[8] 劉羽.基于機群的二維大地電磁Occam反演的并行計算研究[D].武漢:中國地質大學,2006

LIU Y.PC-Cluster based parallel computation research on 2-D magnetotelluric occam inversion[D].Wuhan:China University of Geosciences,2006

[9] 陳露軍.基于MPI的三維井地電磁場并行計算研究與實現(xiàn)[D].成都:成都理工大學,2006

CHEN L J.Research on parallel computation of 3-D borehole-surface EM based on MPI[D].Chengdu:Chengdu University of Technology,2006

[10] SIRIPUNVARAPORN W,EGBERT G.WSINV3D-MT:vertical magnetic field transfer function inversion and parallel implementation[J].Physics of the Earth and Planetary Interiors,2009,173(3/4):317-329

[11] 李焱,胡祥云,吳桂桔,等.基于MPI的二維大地電磁正演的并行計算[J].地震地質,2010,32(3):392-401

LI Y,HU X Y,WU G J,et al.Research of 1-D magnetotelluric parallel computation based on MPI[J].Seismology and Geology,2010,32(3):392-401

[12] CONSTABLE S C,PARKER R L,CONSTABLE C G.Occam’s inversion:a practical algorithm for generating smooth models from EM sounding data[J].Geophysics,1987,52(3):289-300

[13] GROOT-HEDLIN DE C,CONSTABLE S C.Occam’s inversion and the North American Central plains electrical anomaly[J].Journal of Geomagnetismc and Geoelectricity,1993,45(9):985-999

[14] 劉羽.MT Occam 并行反演方案及性能分析[J].武漢理工大學學報,2007,29(12):136-140

LIU Y.Parallel MT occam inversion scheme and its performance analysis[J].Journal of Wuhan University of Technology,2007,29(12):136-140

[15] 馬召貴,趙改善,武港山,等.起伏地表疊前時間偏移的多級并行優(yōu)化技術[J].石油物探,2013,52(3):280-287

MA Z G,ZHAO G S,WU G S,et al.Multi-level parallel optimization technique for prestack time migration from rugged topography[J].Geophysical Prospecting for Petroleum,2013,52(3):280-287

[16] JACOBSEN D A,SENOCAK I.Multi-level parallelism for incompressible flow computations on GPU clusters[J].Parallel Computing,2013,39(1):1-20

[17] 張猛,王華忠,任浩然,等.基于CPU/GPU異構平臺的全波形反演及其實用化分析[J].石油物探,2014,53(4):461-467

ZHANG M,WANG H Z,REN H R,et al.Fall waveform inversion on the CPU/GPU heterogeneous platform and its application analysis[J].Geophysical Prospecting for Petroleum,2014,53(4):461-467

[18] XIAO Y,LIU Y.GPU acceleration for the gaussian elimination in magnetotelluric Occam inversion algorithm[J].Proceedings of the 4thInternational Conference on Computer Engineering and Networks,2015:123-131