魏哲楓,朱成宏,陳業(yè)全

(中國(guó)石油化工股份有限公司石油勘探開(kāi)發(fā)研究院,北京100083)

三維全波形反演高效異構(gòu)并行計(jì)算

魏哲楓,朱成宏,陳業(yè)全

(中國(guó)石油化工股份有限公司石油勘探開(kāi)發(fā)研究院,北京100083)

全波形反演(FWI)每次迭代都需要進(jìn)行若干次地震波正演,計(jì)算量非常大,尤其在三維情況下,提高并行計(jì)算的效率和穩(wěn)健性至關(guān)重要。引入隨機(jī)邊界來(lái)反傳、重建震源波場(chǎng),可充分發(fā)揮GPU的計(jì)算能力,從而實(shí)現(xiàn)反演梯度的高效計(jì)算,相比監(jiān)測(cè)點(diǎn)(checkpoint)和有效邊界技術(shù),大幅減少了數(shù)據(jù)存儲(chǔ)和數(shù)據(jù)交換的開(kāi)銷(xiāo),具有計(jì)算效率高和存儲(chǔ)量小的優(yōu)點(diǎn);開(kāi)發(fā)了作業(yè)池并行作業(yè)管理機(jī)制,與常規(guī)消息傳遞接口(message passing interface,MPI)并行機(jī)制相比,可動(dòng)態(tài)增減節(jié)點(diǎn),具有近似線性的加速比,更適應(yīng)大規(guī)模異構(gòu)并行。采用三維SEG/EAGE推覆體模型進(jìn)行了速度反演測(cè)試,結(jié)果證明該技術(shù)高效且可靠。

全波形反演;隨機(jī)邊界;異構(gòu)并行;作業(yè)池;速度反演

全波形反演(FWI)直接基于波動(dòng)方程,以模型正演與觀測(cè)地震波形一致為依據(jù),可高精度反演速度、密度、波阻抗和Q值等巖性參數(shù)[1-2]。目前,全波形反演精細(xì)速度建模已開(kāi)始應(yīng)用于油氣勘探地震資料成像中[3-5],但大規(guī)模工業(yè)化應(yīng)用仍面臨諸多挑戰(zhàn)。

局部極小值是FWI這類(lèi)非線性反演的固有難題,為此,學(xué)者們提出了許多提高FWI目標(biāo)函數(shù)線性度的方法[6-9]。海量計(jì)算是FWI的另一個(gè)難題,尤其是3D FWI,其計(jì)算開(kāi)銷(xiāo)是常規(guī)逆時(shí)偏移的數(shù)倍乃至數(shù)十倍,因此基于GPU集群的高性能計(jì)算技術(shù)受到廣泛的重視[10-11]。

3D FWI高性能計(jì)算方案的核心是實(shí)現(xiàn)單炮反演梯度的高效計(jì)算和多炮作業(yè)間的穩(wěn)健并行管理。FWI的梯度計(jì)算需要震源波場(chǎng)和檢波點(diǎn)波場(chǎng)同時(shí)刻相關(guān),類(lèi)似于完成一次逆時(shí)偏移(RTM)[12-13],因此可借鑒RTM的不同實(shí)現(xiàn)方法,包括:完全保存震源波場(chǎng)的完全存儲(chǔ)法;優(yōu)選有限波場(chǎng)時(shí)間切片用于反傳重建震源波場(chǎng)的監(jiān)測(cè)點(diǎn)方法[14-15];保存所有時(shí)刻邊界波場(chǎng)值以反傳重建震源波場(chǎng)的有效邊界法[16];引入隨機(jī)邊界,保存最后兩個(gè)時(shí)刻的隨機(jī)波場(chǎng)以實(shí)現(xiàn)震源波場(chǎng)重建的隨機(jī)邊界法[17]。一般來(lái)說(shuō),GPU存儲(chǔ)空間相對(duì)有限,難以滿足完全存儲(chǔ)法和監(jiān)測(cè)點(diǎn)方法對(duì)波場(chǎng)存儲(chǔ)的要求,如果借助主機(jī)端內(nèi)存或硬盤(pán)來(lái)保存,則會(huì)增加GPU端與CPU端的數(shù)據(jù)交換,極大地犧牲GPU計(jì)算時(shí)間。有效邊界方法在進(jìn)行邊界賦值時(shí)難以滿足GPU合并訪存要求,會(huì)帶來(lái)高延遲,降低GPU計(jì)算效率。隨機(jī)邊界法只需要保存最后兩個(gè)時(shí)刻的隨機(jī)波場(chǎng),存儲(chǔ)量小(當(dāng)前的主流GPU顯存已能滿足存儲(chǔ)需要),且隨機(jī)波場(chǎng)的載入很容易滿足GPU合并訪存要求,能充分發(fā)揮GPU的高效計(jì)算能力。

FWI多炮并行作業(yè)管理的常用方案是消息傳遞接口(message passing interface,MPI)[10],利用MPI提供的消息通信接口,將并行集群中參與計(jì)算的諸多節(jié)點(diǎn)納入統(tǒng)一的通信域中,節(jié)點(diǎn)間通過(guò)消息傳遞實(shí)現(xiàn)數(shù)據(jù)交換、同步、互斥等協(xié)同工作?；贛PI的節(jié)點(diǎn)間并行實(shí)現(xiàn)在編程上更簡(jiǎn)潔直觀,有助于減輕開(kāi)發(fā)人員的工作強(qiáng)度,但是其缺點(diǎn)是[18]:對(duì)于大規(guī)模節(jié)點(diǎn)共同參與的并行作業(yè)一旦某一個(gè)節(jié)點(diǎn)宕機(jī)則會(huì)導(dǎo)致整個(gè)作業(yè)終止。此外,由于節(jié)點(diǎn)間通信會(huì)導(dǎo)致網(wǎng)絡(luò)阻塞,并行節(jié)點(diǎn)規(guī)模越大,MPI并行加速比越低。

本文首先簡(jiǎn)單回顧時(shí)間域FWI的基本理論;然后介紹為提高梯度計(jì)算效率,采用隨機(jī)邊界實(shí)現(xiàn)震源波場(chǎng)高效重建的方法;接著重點(diǎn)介紹為克服MPI多節(jié)點(diǎn)并行作業(yè)管理的缺點(diǎn),基于共享存儲(chǔ)空間和文件鎖技術(shù)的“作業(yè)池”機(jī)制的設(shè)計(jì)與實(shí)現(xiàn);最后用三維SEG/EAGE推覆體模型進(jìn)行數(shù)值測(cè)試,驗(yàn)證所提出的方法技術(shù)的正確性。

1 方法原理

全波形反演每次迭代都需要若干次地震波場(chǎng)模擬,為便于處理介質(zhì)的非均勻性,使用如下聲波方程[19]:

(1)

式中:R3表示三維實(shí)空間;p(x,t)為波場(chǎng)值;ρ(x)和v(x)分別為地下介質(zhì)密度和速度;xs為炮點(diǎn)坐標(biāo);s(xs,t)為震源函數(shù)。定義p(xr,t;xs)為給定模型的數(shù)值模擬波場(chǎng)值,d(xr,t;xs)為實(shí)際地震記錄,波場(chǎng)殘差δp(xr,t;xs)=p(xr,t;xs)-d(xr,t;xs),其中,xr為檢波點(diǎn)坐標(biāo)。則反演的目標(biāo)函數(shù)可表示為:

(2)

全波形反演就是尋求目標(biāo)函數(shù)S(v)的極小值,滿足目標(biāo)函數(shù)極小值的速度即為反演的最終速度。由Born近似理論及伴隨方法[4-5]可推導(dǎo)出目標(biāo)函數(shù)S(v)對(duì)速度v的導(dǎo)數(shù)為:

(3)

其中,λ(x,t;xr)為波場(chǎng)殘差反傳波場(chǎng),滿足:

(4)

此反傳波場(chǎng)方程稱為伴隨方程,它以波場(chǎng)殘差作為震源。

1.1 反演梯度的高效計(jì)算

為了在GPU/CPU異構(gòu)的高性能計(jì)算平臺(tái)上實(shí)現(xiàn)3D FWI梯度計(jì)算并得到很高的計(jì)算效率,我們將CLAPP提出的隨機(jī)邊界條件[17]與GPU加速相結(jié)合,高效重建震源波場(chǎng)。隨機(jī)邊界能夠保存波場(chǎng)能量,同時(shí)又打破波前面的相干性,利用該特性可以重建震源波場(chǎng),并且進(jìn)一步將其與檢波點(diǎn)波場(chǎng)同時(shí)刻相關(guān),得到反演梯度。與有效邊界存儲(chǔ)法和監(jiān)測(cè)點(diǎn)方法相比,這種方法存儲(chǔ)量小、計(jì)算效率高。

在邊界區(qū)域設(shè)置隨機(jī)變化的速度場(chǎng)構(gòu)成隨機(jī)邊界條件:

(5)

式中:α為與速度有關(guān)的量,與速度有相同的量綱;R代表隨機(jī)邊界區(qū)域的厚度;r代表隨機(jī)邊界區(qū)域中的點(diǎn)距邊界的距離;Ftaper∈(0,1)為一平滑窗函數(shù);Ran(Sseed)代表由種子數(shù)Sseed產(chǎn)生的隨機(jī)數(shù),Ran (Sseed)∈(0,1);xb代表隨機(jī)邊界中的點(diǎn);v(xb)代表隨機(jī)邊界中的點(diǎn)xb處的隨機(jī)速度值;v0(xb)代表隨機(jī)邊界上r=0處的速度值,由它產(chǎn)生隨機(jī)邊界中各點(diǎn)的隨機(jī)速度值。通過(guò)調(diào)節(jié)α和R可以得到合適的隨機(jī)邊界。

以速度為2500m/s的均勻介質(zhì)為例,如圖1所示,在模型邊界處設(shè)置50個(gè)網(wǎng)格層的隨機(jī)邊界(網(wǎng)格間隔為20m),R=50,v0=2500m/s,α=1000m/s。設(shè)置震源在介質(zhì)中心位置,先正傳波場(chǎng)到隨機(jī)邊界上,波場(chǎng)會(huì)被邊界處的隨機(jī)介質(zhì)散射打亂,再將打散的地震波場(chǎng)反傳,即可重建任一過(guò)去時(shí)刻的震源波場(chǎng)。圖2給出了某一時(shí)刻的反傳重建波場(chǎng)快照,可見(jiàn)波場(chǎng)呈現(xiàn)散亂的隨機(jī)狀,說(shuō)明邊界處介質(zhì)的隨機(jī)性設(shè)置得較適當(dāng)。

隨機(jī)邊界反傳重建波場(chǎng)在效率和存儲(chǔ)方面具有明顯優(yōu)勢(shì),但也有不足之處:重建波場(chǎng)中的低頻成分會(huì)有損失,不過(guò)仍能滿足成像精度的需要,前人的工作證明了這一點(diǎn)[17]。

1.2 大規(guī)模并行作業(yè)管理機(jī)制

全波形反演的每次迭代計(jì)算包含計(jì)算梯度、線性搜索、更新速度模型等步驟,而計(jì)算梯度、線性搜索等過(guò)程的實(shí)現(xiàn)(圖3),需要進(jìn)行大量炮點(diǎn)的地震波場(chǎng)模擬[20]。3D FWI需要大型計(jì)算機(jī)集群并行運(yùn)算,如何實(shí)現(xiàn)靈活、穩(wěn)健和高效的全波形反演的大規(guī)模并行作業(yè)管理,是實(shí)現(xiàn)產(chǎn)業(yè)化全波形反演技術(shù)的一大難點(diǎn)。

圖1 隨機(jī)邊界

圖2 利用隨機(jī)邊界重建震源波場(chǎng)

目前已公開(kāi)的全波形反演并行作業(yè)管理方案,都是采用MPI技術(shù)實(shí)現(xiàn)(圖4),參與計(jì)算的節(jié)點(diǎn)之間通過(guò)MPI互相關(guān)聯(lián)在一個(gè)通信域中,節(jié)點(diǎn)之間通過(guò)網(wǎng)絡(luò)通訊交換數(shù)據(jù)。MPI編程有兩個(gè)固有缺點(diǎn):①不支持計(jì)算節(jié)點(diǎn)的動(dòng)態(tài)增減;②任何一個(gè)計(jì)算節(jié)點(diǎn)如出現(xiàn)故障,則整個(gè)通信域中的節(jié)點(diǎn)都會(huì)停止作業(yè)?；贛PI的FWI并行計(jì)算實(shí)現(xiàn)方法通常有3種(為便于說(shuō)明,這里以20炮為例,數(shù)字1～20分別表示1～20炮地震數(shù)據(jù),A,B,C分別表示參與計(jì)算的GPU/CPU異構(gòu)集群節(jié)點(diǎn),假定A和B的計(jì)算效率相同,C的計(jì)算效率是A或B的兩倍,以下討論的并行作業(yè)管理是關(guān)于節(jié)點(diǎn)間或進(jìn)程間的,因此,對(duì)于單純的CPU集群同樣適用)。第1種是基于單炮區(qū)域分解的MPI并行實(shí)現(xiàn)方法。如圖5所示(圖中編號(hào)1～20表示炮號(hào),大箭頭僅表示計(jì)算順序,無(wú)輸入、輸出的依賴關(guān)系),該方法依次計(jì)算炮1,2,3,…,20,并在計(jì)算每一炮的過(guò)程中,把該炮計(jì)算區(qū)域按參與計(jì)算的節(jié)點(diǎn)數(shù)(或進(jìn)程數(shù))分解,所有節(jié)點(diǎn)都將參與每一炮的計(jì)算,在計(jì)算過(guò)程中,所有計(jì)算節(jié)點(diǎn)(或進(jìn)程)需要在每一個(gè)計(jì)算時(shí)刻點(diǎn)同步,并相互交換邊界重疊區(qū)域的計(jì)算結(jié)果。其優(yōu)點(diǎn)是:多節(jié)點(diǎn)協(xié)作能夠克服單個(gè)計(jì)算節(jié)點(diǎn)內(nèi)存不足的問(wèn)題。其缺點(diǎn)是:①計(jì)算效率高的節(jié)點(diǎn)(或進(jìn)程)總是在每一個(gè)時(shí)刻點(diǎn)等待計(jì)算效率低的節(jié)點(diǎn),導(dǎo)致整體并行計(jì)算效率下降;②節(jié)點(diǎn)之間需要在每個(gè)時(shí)刻交換數(shù)據(jù),因而該方法受計(jì)算節(jié)點(diǎn)之間網(wǎng)絡(luò)帶寬限制。第2種是基于均分炮數(shù)的MPI并行實(shí)現(xiàn)方法。如圖6 所示,20炮被盡可能均勻地分配到每一個(gè)計(jì)算節(jié)點(diǎn)上,每一個(gè)計(jì)算節(jié)點(diǎn)負(fù)責(zé)一部分炮的計(jì)算。其優(yōu)點(diǎn)是:計(jì)算節(jié)點(diǎn)之間無(wú)需交換數(shù)據(jù),不受網(wǎng)絡(luò)帶寬的限制。其缺點(diǎn)是:計(jì)算效率高的節(jié)點(diǎn)要等待計(jì)算效率最低的節(jié)點(diǎn),導(dǎo)致整體并行計(jì)算效率下降。第3種方法是基于主從編程模型,存在一個(gè)主節(jié)點(diǎn)負(fù)責(zé)并行作業(yè)的調(diào)度和管理,從節(jié)點(diǎn)通過(guò)與主節(jié)點(diǎn)通信來(lái)獲取計(jì)算任務(wù)。其優(yōu)點(diǎn)是:可以充分利用計(jì)算節(jié)點(diǎn)的計(jì)算能力,保證能者多勞。其缺點(diǎn)是:海量節(jié)點(diǎn)的大規(guī)模并行計(jì)算會(huì)出現(xiàn)通信擁堵,因?yàn)樘嗟膹墓?jié)點(diǎn)會(huì)在同一時(shí)刻要求與主節(jié)點(diǎn)通信,會(huì)對(duì)網(wǎng)絡(luò)帶來(lái)極大壓力,使得通信成本太高,計(jì)算效率無(wú)法隨著節(jié)點(diǎn)數(shù)的增加而線性提高。除了3種基本實(shí)現(xiàn)方法,還可以將3種方法相互組合形成混合式的并行作業(yè)管理方案,在此不再一一贅述。

圖3 全波形反演流程

圖4 全波形反演MPI并行計(jì)算實(shí)現(xiàn)過(guò)程

綜上所述,在FWI技術(shù)實(shí)現(xiàn)中,采用基于MPI的并行作業(yè)管理方案存在諸多不足:①不支持計(jì)算節(jié)點(diǎn)的動(dòng)態(tài)增加和減少,缺乏靈活性,難以適應(yīng)實(shí)際生產(chǎn)作業(yè)中的需求變動(dòng);②穩(wěn)定性差,采用MPI技術(shù)時(shí),即使節(jié)點(diǎn)之間不發(fā)生數(shù)據(jù)交換,也相互關(guān)聯(lián),一旦某個(gè)節(jié)點(diǎn)出現(xiàn)硬件故障,整個(gè)計(jì)算立即失敗;③對(duì)異構(gòu)設(shè)備適應(yīng)性差,異構(gòu)的計(jì)算機(jī)集群各節(jié)點(diǎn)間效率不一致,采用MPI并行實(shí)現(xiàn)會(huì)導(dǎo)致整體計(jì)算效率將被最慢的節(jié)點(diǎn)拖累;④硬件成本高,MPI并行實(shí)現(xiàn)對(duì)計(jì)算機(jī)集群節(jié)點(diǎn)之間的網(wǎng)絡(luò)速度和穩(wěn)定性要求很高,迫使人們采用Infiniband以及其它較昂貴的設(shè)備來(lái)提高網(wǎng)絡(luò)通信效率和穩(wěn)定性;⑤大規(guī)模并行計(jì)算效率低下,MPI技術(shù)的并行效率隨計(jì)算機(jī)節(jié)點(diǎn)數(shù)增加而降低(尤其當(dāng)節(jié)點(diǎn)數(shù)達(dá)到幾十個(gè)以上時(shí)),無(wú)法達(dá)到近似線性的加速比,MPI對(duì)計(jì)算效率提升的瓶頸作用就有體現(xiàn)。

圖5 基于單炮區(qū)域分解的MPI并行實(shí)現(xiàn)方法

圖6 基于均分炮數(shù)的MPI并行實(shí)現(xiàn)方法

為克服全波形反演技術(shù)采用MPI并行實(shí)現(xiàn)中存在的以上不足,我們研發(fā)了基于“作業(yè)池”的并行作業(yè)管理機(jī)制,具有如下技術(shù)優(yōu)勢(shì):并行計(jì)算節(jié)點(diǎn)的動(dòng)態(tài)增減,具備硬件容錯(cuò)能力,支持異構(gòu)設(shè)備,對(duì)網(wǎng)絡(luò)等硬件設(shè)施要求低,并行計(jì)算效率隨計(jì)算節(jié)點(diǎn)近似線性增加(這里“節(jié)點(diǎn)”可以是物理上的計(jì)算節(jié)點(diǎn),也可以是并行計(jì)算中的進(jìn)程,下同)。該方法對(duì)每一個(gè)全波形反演項(xiàng)目,首先根據(jù)項(xiàng)目的具體情況產(chǎn)生一個(gè)作業(yè)狀態(tài)文件,如圖7所示,用來(lái)記錄每一次并行計(jì)算過(guò)程中每一炮的計(jì)算狀態(tài)。這些計(jì)算狀態(tài)分為:已經(jīng)完成(Done)、沒(méi)有被計(jì)算(Untouched)、有一個(gè)節(jié)點(diǎn)在計(jì)算(1stTry)、有兩個(gè)節(jié)點(diǎn)在計(jì)算(2ndTry)、有3個(gè)節(jié)點(diǎn)在計(jì)算(3rdTry)……。在每一次并行計(jì)算過(guò)程中,采用自動(dòng)作業(yè)控制方法來(lái)實(shí)現(xiàn)“按需分配”計(jì)算機(jī)資源,如圖8所示,計(jì)算節(jié)點(diǎn)A,B,C在每次并行計(jì)算時(shí),會(huì)依次查看所有的炮(1～20),如果本次并行中某炮還是Untouched狀態(tài),該節(jié)點(diǎn)就啟動(dòng)對(duì)該炮的計(jì)算,并把本次并行中該炮的狀態(tài)更改為1stTry狀態(tài),并且一旦計(jì)算完成,就把該炮狀態(tài)更改為Done狀態(tài)。如果一個(gè)計(jì)算節(jié)點(diǎn)發(fā)現(xiàn)本次并行計(jì)算中沒(méi)有任何一炮在Untouched狀態(tài),就依次查看是否有1stTry狀態(tài)的炮,如果有,該節(jié)點(diǎn)立即啟動(dòng)對(duì)該炮的計(jì)算,并把本次并行中該炮的狀態(tài)更改為2ndTry狀態(tài)。依次類(lèi)推,直到所有炮都被計(jì)算完而且狀態(tài)都為Done狀態(tài),本次并行完成。在具體編程實(shí)現(xiàn)上,為使不同節(jié)點(diǎn)對(duì)作業(yè)狀態(tài)文件的操作不出現(xiàn)競(jìng)態(tài)條件問(wèn)題(race condition),我們需要對(duì)作業(yè)狀態(tài)文件進(jìn)行加鎖處理(例如,利用fcntl來(lái)鎖文件)?；谧詣?dòng)作業(yè)控制的并行實(shí)現(xiàn)方法能夠充分實(shí)現(xiàn)異構(gòu)計(jì)算機(jī)集群的高效并行,并且任何節(jié)點(diǎn)如果出現(xiàn)硬件故障(如死機(jī)),都不會(huì)導(dǎo)致整個(gè)全波形反演項(xiàng)目中斷。

每個(gè)計(jì)算節(jié)點(diǎn)都是獨(dú)立的主節(jié)點(diǎn)(主進(jìn)程),共用一個(gè)共享盤(pán),節(jié)點(diǎn)和節(jié)點(diǎn)之間沒(méi)有直接通信,節(jié)點(diǎn)與共享盤(pán)之間只有在控制點(diǎn)上有數(shù)據(jù)交換。在全波形反演計(jì)算過(guò)程中,每個(gè)計(jì)算節(jié)點(diǎn)上都獨(dú)立啟動(dòng)一個(gè)作業(yè)(主進(jìn)程),這些主進(jìn)程會(huì)通過(guò)競(jìng)爭(zhēng)的方式讀寫(xiě)位于共享盤(pán)中的狀態(tài)文件和共享數(shù)據(jù)體,如圖9右邊所示。所有主進(jìn)程都能通過(guò)作業(yè)狀態(tài)文件找到參與全波形反演計(jì)算的切入點(diǎn),修改作業(yè)狀態(tài)文件。這種分布式并行計(jì)算節(jié)點(diǎn)之間沒(méi)有緊密的依存關(guān)系,因此在實(shí)際應(yīng)用中具有如下優(yōu)點(diǎn):①支持硬件設(shè)備的動(dòng)態(tài)增加和減少,參與計(jì)算的節(jié)點(diǎn)可以動(dòng)態(tài)增加和減少,每個(gè)節(jié)點(diǎn)/進(jìn)程都可以通過(guò)共享狀態(tài)文件/數(shù)據(jù)來(lái)調(diào)節(jié)工作切入點(diǎn);②穩(wěn)定性好,由于計(jì)算機(jī)節(jié)點(diǎn)無(wú)相互關(guān)聯(lián),在長(zhǎng)時(shí)間的全波形反演計(jì)算中,即使部分節(jié)點(diǎn)出現(xiàn)故障中斷,整個(gè)計(jì)算作業(yè)仍會(huì)繼續(xù);③支持異構(gòu)設(shè)備,不同計(jì)算效率的節(jié)點(diǎn)能夠根據(jù)狀態(tài)文件自己調(diào)整計(jì)算量,不會(huì)干預(yù)其它節(jié)點(diǎn)的計(jì)算;④硬件成本低,并行計(jì)算對(duì)系統(tǒng)的穩(wěn)定性、網(wǎng)絡(luò)通信等要求比較低,因此可以使用低成本的硬件設(shè)備;⑤并行計(jì)算效率不受節(jié)點(diǎn)數(shù)目影響,由于沒(méi)有節(jié)點(diǎn)間的網(wǎng)絡(luò)通訊瓶頸,增加節(jié)點(diǎn)不會(huì)帶來(lái)額外的計(jì)算效率損失。

圖7 并行作業(yè)狀態(tài)控制文件

圖8 基于自動(dòng)作業(yè)控制的并行實(shí)現(xiàn)方法

圖9 全波形反演分布式并行計(jì)算實(shí)現(xiàn)方法

2 應(yīng)用算例

三維SEG/EAGE推覆體模型是檢驗(yàn)三維全波形反演方法正確性的通用模型,模型大小為20.00km×20.00km×4.65km,速度范圍為2200～6000m/s(圖10)。該模型模仿野外復(fù)雜地質(zhì)構(gòu)造,主要由受構(gòu)造運(yùn)動(dòng)影響拉張分裂的基巖和推覆于其上的沉積巖序列兩部分構(gòu)成。其頂界面為一不整合面。在此不整合面之上,發(fā)育有12層連續(xù)沉積層序,每層厚度設(shè)計(jì)尺度與野外實(shí)際觀測(cè)相同,沉積序列被其底部發(fā)育的巖體分開(kāi)。推覆體模型的巖相及地層學(xué)特征模擬了區(qū)域張應(yīng)力作用下的海相與構(gòu)造混合沉積分布于某些沉積層位中的河道及裂縫。在此區(qū)域中,所有封閉的圈閉均含有油氣。

推覆體模型構(gòu)造復(fù)雜,縱向?qū)游欢嘧儾⒕哂袕?qiáng)烈的橫向速度變化,可以用來(lái)檢驗(yàn)全波形反演方法的適應(yīng)性和其對(duì)復(fù)雜構(gòu)造的成像能力。斷層的平面分布等構(gòu)造現(xiàn)象可以考察全波形反演對(duì)構(gòu)造細(xì)節(jié)刻畫(huà)的能力。尤為重要的是,中深層較小的含油氣圈閉,可以有效地反映全波形反演的高分辨能力以及儲(chǔ)層刻畫(huà)能力。

采用圖11所示的平滑模型作為反演的初始速度模型,設(shè)置10000炮用于反演測(cè)試,沿inline和xline等間隔分別布設(shè)100炮,雙邊放炮,最大偏移距8000m,道間距40m,反演前先對(duì)數(shù)據(jù)進(jìn)行帶通濾波,濾波后主頻為8Hz。計(jì)算平臺(tái)是20個(gè)配有GPU M2090節(jié)點(diǎn)和28個(gè)配有GPU K20節(jié)點(diǎn)組成的GPU/CPU異構(gòu)高性能計(jì)算機(jī)集群,且集群節(jié)點(diǎn)的計(jì)算能力不同(主要體現(xiàn)在GPU卡M2090和K20之間的計(jì)算效率有差異),集群節(jié)點(diǎn)的資源配置見(jiàn)表1和表2,反演的計(jì)算規(guī)模和效率見(jiàn)表3。作業(yè)池并行機(jī)制實(shí)現(xiàn)了節(jié)點(diǎn)間的解耦,整個(gè)作業(yè)調(diào)度過(guò)程采用了競(jìng)爭(zhēng)機(jī)制,所有計(jì)算節(jié)點(diǎn)到“作業(yè)池”中搶任務(wù),直到所有任務(wù)完成為止,相比傳統(tǒng)MPI并行機(jī)制既減少了通信開(kāi)銷(xiāo),又不會(huì)被效率低的節(jié)點(diǎn)所拖累。另外,為了測(cè)試作業(yè)池并行機(jī)制的穩(wěn)健性,我們?cè)诜囱葸^(guò)程中會(huì)每隔若干小時(shí)人為停止某些進(jìn)程,并相應(yīng)地加入一些新節(jié)點(diǎn),以模擬硬件故障情況下程序并行機(jī)制的穩(wěn)健性。經(jīng)過(guò)40次反演,程序運(yùn)行穩(wěn)健,由不同配置節(jié)點(diǎn)組成的集群在動(dòng)態(tài)增減節(jié)點(diǎn)的情況下順利完成反演,運(yùn)算中各節(jié)點(diǎn)根據(jù)各自的計(jì)算能力自動(dòng)實(shí)現(xiàn)均衡負(fù)載,總用時(shí)約67h。圖12a,圖12b和圖12c分別為深度在1450m處,inline和xline均位于11200m處的初始速度模型、迭代40次的反演速度模型和真實(shí)速度模型切片。與初始模型相比,反演模型河道識(shí)別清晰,小的速度異常體也得到呈現(xiàn),斷層斷面刻畫(huà)清楚。圖13是迭代40次對(duì)應(yīng)的反演梯度,其對(duì)應(yīng)了此次迭代速度模型的更新量。可見(jiàn),基于作業(yè)池并行機(jī)制的FWI技術(shù)實(shí)用、可靠,并行節(jié)點(diǎn)間的解耦使得節(jié)點(diǎn)間互不干擾,既靈活又穩(wěn)健,競(jìng)爭(zhēng)式的作業(yè)調(diào)度提高了并行效率。

圖10 三維SEG/EAGE推覆體模型

圖11 三維SEG/EAGE推覆體初始平滑模型

操作系統(tǒng)RedHatEnterpriseLinuxServerrelease5.8CUDADriverVersion/RuntimeVersion6.5/5.5CPU型號(hào)Intel(R)Xeon(R)CPUE5-26600@2.20GHzCPU核數(shù)16內(nèi)存/GB62每個(gè)節(jié)點(diǎn)上的GPU卡數(shù)2塊

表2 GPU設(shè)備的參數(shù)

圖12 深度在1450m處,inline和xline均位于11200m處的速度模型切片a 初始速度模型; b 迭代40次的反演速度模型; c 真實(shí)速度模型

動(dòng)用的計(jì)算資源20個(gè)配備M2090的節(jié)點(diǎn);28個(gè)配備K20的節(jié)點(diǎn)總炮數(shù)/參與反演炮數(shù)10000/5000單炮網(wǎng)格點(diǎn)數(shù)18210801(301×301×201)單炮時(shí)間步數(shù)2501單炮正演平均時(shí)間/s15.18(M2090);8.71(K20)單次迭代平均用時(shí)/h1.68

圖13 迭代40次的反演梯度

3 結(jié)論

本文分析了三維全波形反演中面臨的梯度計(jì)算效率及并行穩(wěn)定性難題,提出將GPU加速技術(shù)與隨機(jī)邊界相結(jié)合以構(gòu)建梯度計(jì)算引擎,彌補(bǔ)了GPU在數(shù)據(jù)存儲(chǔ)和交換方面的不足,充分發(fā)揮了GPU加速技術(shù)的優(yōu)勢(shì);研發(fā)的作業(yè)池機(jī)制,利用共享的作業(yè)狀態(tài)文件控制并行作業(yè)的進(jìn)程,避免了計(jì)算節(jié)點(diǎn)間的直接通信和相互影響,提高了大規(guī)模并行計(jì)算的穩(wěn)健性。

三維SEG/EAGE推覆體模型算例驗(yàn)證了本文并行計(jì)算技術(shù)的高效性和穩(wěn)健性,作業(yè)池機(jī)制帶來(lái)的計(jì)算節(jié)點(diǎn)動(dòng)態(tài)增減功能,增強(qiáng)了FWI對(duì)不同硬件條件和計(jì)算資源變動(dòng)的適應(yīng)性,有助于全波形反演技術(shù)的大規(guī)模工業(yè)化應(yīng)用。

[1] LAILLY P.The seismic inverse problem as a sequence of before stack migrations[J].Expanded Abstracts of conference on inverse scattering,theory and application,society for industrial and applied mathematics,1983:206-220

[2] TARANTOLA A.Inversion of seismic-reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

[3] KAPOOR S,VIGH D,WIARDA E,et al.Full waveform inversion around the world[J].Expanded Abstracts of 75thEAGE Conference&Exhibition incorporating SPE,2013:We-11-03

[4] JOSEPH M,TIMOTHY K,THIERRY T,et al.3D acoustic waveform inversion of land data:a case study from Saudi Arabia[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1-5

[5] PLESSIX R,BAETEN G,DEMAAG J W,et al.Full waveform inversion and distance separated simultaneous sweeping:a study with a land seismic data set[J].Geophysical Prospecting,2012,60(4):733-747

[6] LUO Y,SCHUSTER G T.Wave-equation traveltime inversion[J].Geophysics,1991,56(5):645-653

[7] SHIN C,CHA Y H.Waveform inversion in the Laplace domain[J].Geophysical Journal International,2008,173(3):922-931

[8] XU S,WANG D L,CHEN F,et al.Full waveform inversion for reflected seismic data[J].Expanded Abstracts of 74thEAGE Conference & Exhibition incorporating SPE EUROPEC,2012:W024

[10] YOUNGSEO K,CHANGSOO S,HENRI C,et al.An algorithm for 3D acoustic time-Laplace-Fourier-domain hybrid full waveform inversion[J].Geophysics,2013,78(4):R151-R166

[11] YANG P L,GAO J H,WANG B L.A graphics processing unit implementation of time-domain full-waveform inversion[J].Geophysics,2015,80(3):F31-F39

[12] CHAVENT G,PLESSIX R E.An optimal true amplitude least-squares prestack depth-migration operator[J].Geophysics,1999,64(2):508-515

[13] SHIN C,MIN D J,YANG D,et al.Evaluation of poststack migration in terms of virtual source and partial de-rivative wavefields[J].Journal of Seismic Exploration,2003,12(1):17-37

[14] SYMES W.Reverse time migration with optimal check pointing[J].Geophysics,2007,72(5):SM213-SM221

[15] ANDERSON J E,TAN L J,WANG D.Time-reversal checkpointing methods for RTM and FWI[J].Geophysics,2012,77(4):S93-S103

[16] GAO J H,YANG P L,WANG B L.Using the effective boundary saving strategy in GPU-based RTM programming[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:4014-4021

[17] CLAPP R G.Reverse time migration with random boundaries[J].Expanded Abstracts of 79thAnnual Internat SEG Mtg,2009:2809-2813

[18] SUH Y,YEH A,WANG B,et al.Cluster programming for reverse time migration[J].The Leading Edge,2010,29(1):94-97

[19] 徐義,張劍鋒．地震波數(shù)值模擬的非規(guī)則網(wǎng)格PML吸收邊界[J]．地球物理學(xué)報(bào),2008,51(5):1520-1526 XU Y,ZHANG J F.An irregular-grid perfectly matched layer absorbing boundary for seismic wave modeling[J].Chinese Journal of Geophysics,2008,51(5):1520-1526

[20] 胡光輝,賈春梅,夏洪瑞,等.三維聲波全波形反演的實(shí)現(xiàn)與驗(yàn)證[J].石油物探,2013,52(4):417-425 HU G H,JIA C M,XIA H R,et al.Implementation and validation of 3D acoustic full waveform inversion[J].Geophysical Prospecting for Petroleum,2013,52(4):417-425

(編輯：陳 杰)

Efficient heterogeneous parallel computing of 3D full waveform inversion

WEI Zhefeng,ZHU Chenghong,CHEN Yequan

(PetroleumExplorationandProductionResearchInstitute,SINOPEC,Beijing100083,China)

Each iteration of full waveform inversion (FWI) requires several seismic forward modeling and takes large computational effort especially in 3D cases.It is key to improve the efficiency and robustness of parallel computing for FWI.To maximize the ability of GPU,we introduce random boundary to back propagate source wavefield for reconstruction to efficiently compute the gradient inversion.Compared to checkpoints and effective boundary technology,it can significantly reduce the cost of data storage and exchange.To adapt to massive heterogeneous parallel clusters,we develop a parallel job management mechanism named job pool.Compared to the conventional message passing interface (MPI) mechanism,it can add in or remove out computing nodes on the fly and the speedup ratio is approximately linear.Numerical example of 3D Overthrust model for velocity inversion showed high efficiency and reliability of the method.

full waveform inversion,random boundary,heterogeneous parallel,job pool,velocity modeling

2016-10-08;改回日期：2016-12-01。

魏哲楓(1985—),男,博士,主要從事高精度地震速度建模與成像研究。

國(guó)家科技重大專(zhuān)項(xiàng)(2011ZX05031-001-04)資助。

P631

A

1000-1441(2017)01-0089-10

10.3969/j.issn.1000-1441.2017.01.011

This research is financially supported by the National Science and Technology Major Project of China (Grant No.2011ZX05031-001-04 ).