高新成 (東北石油大學(xué)現(xiàn)代教育技術(shù)中心，黑龍江 大慶 163318)石穎 (東北石油大學(xué)地球科學(xué)學(xué)院，黑龍江 大慶 163318)

?

基于CPU/GPU異構(gòu)平臺(tái)的疊前逆時(shí)偏移成像系統(tǒng)

高新成 (東北石油大學(xué)現(xiàn)代教育技術(shù)中心，黑龍江 大慶 163318)石穎 (東北石油大學(xué)地球科學(xué)學(xué)院，黑龍江 大慶 163318)

針對(duì)目前地震資料偏移成像過(guò)程中代碼復(fù)雜分散、可移植性差和可視化操作弱等問(wèn)題，結(jié)合疊前逆時(shí)偏移業(yè)務(wù)需求，研究了CPU/GPU協(xié)同并行計(jì)算技術(shù)優(yōu)化疊前逆時(shí)偏移成像算法，設(shè)計(jì)了逆時(shí)偏移系統(tǒng)流程與功能，開(kāi)發(fā)了一套基于CPU/GPU異構(gòu)平臺(tái)的疊前逆時(shí)偏移成像系統(tǒng)。利用自設(shè)模型和傳統(tǒng)模型數(shù)據(jù)對(duì)系統(tǒng)進(jìn)行測(cè)試，測(cè)試結(jié)果表明系統(tǒng)操作方便可行，運(yùn)行效果良好，能夠有效的對(duì)復(fù)雜構(gòu)造地質(zhì)體精確成像。

逆時(shí)偏移；CPU/GPU并行加速；PML邊界條件；QT

隨著地震資料偏移成像應(yīng)用推廣，在逆時(shí)偏移成像處理中積累了很多的應(yīng)用程序和方法，如何將這些程序和方法集成整合起來(lái)，實(shí)現(xiàn)偏移成像的可視化操作和高效管理成為當(dāng)前研究的重要課題[1]。目前疊前逆時(shí)偏移是獲得復(fù)雜構(gòu)造精確成像最有效的方法[2]。為此，筆者結(jié)合疊前逆時(shí)偏移業(yè)務(wù)需求研究地震波場(chǎng)正演模擬、成像條件、邊界吸收條件和噪音壓制算法等，形成了一套疊前逆時(shí)偏移成像方法；并研究了CPU/GPU異構(gòu)平臺(tái)并行加速技術(shù)，將其應(yīng)用于疊前逆時(shí)偏移成像方法中地震波場(chǎng)正演計(jì)算、成像條件和波場(chǎng)反傳計(jì)算等計(jì)算環(huán)節(jié)，開(kāi)發(fā)了一套基于CPU/GPU異構(gòu)平臺(tái)的疊前逆時(shí)偏移成像系統(tǒng)，為油田開(kāi)發(fā)資料的解釋提供技術(shù)支持。

1 CPU/GPU異構(gòu)平臺(tái)的逆時(shí)偏移流程

1.1 疊前逆時(shí)偏移成像算法

筆者采用逆時(shí)偏移算法，在波場(chǎng)正演過(guò)程中采用隨機(jī)邊界條件和PML完全吸收邊界條件相結(jié)合的2種邊界處理方式，成像過(guò)程中采用互相關(guān)成像條件，成像后采用相應(yīng)濾波法進(jìn)行去噪[3～7]。逆時(shí)偏移算法實(shí)現(xiàn)流程如圖1所示。算法實(shí)現(xiàn)步驟如下：

步1 使用波動(dòng)方程有限差分法沿時(shí)間正向延拓波場(chǎng)計(jì)算到Tmax時(shí)刻，保存最后2個(gè)時(shí)刻波場(chǎng)值。

步2 采用PML邊界條件，使用波動(dòng)方程有限差分法沿時(shí)間正向延拓波場(chǎng)計(jì)算到Tmax時(shí)刻，保存每個(gè)時(shí)刻的檢波點(diǎn)數(shù)據(jù)。

步3 讀取最后2個(gè)時(shí)刻波場(chǎng)值，使用波動(dòng)方程有限差分法沿時(shí)間反向延拓計(jì)算到0時(shí)刻，計(jì)算每個(gè)時(shí)刻波場(chǎng)值。

步4 讀取檢波點(diǎn)數(shù)據(jù)，使用波動(dòng)方程有限差分法沿時(shí)間反向延拓計(jì)算到0時(shí)刻，計(jì)算每個(gè)時(shí)刻波場(chǎng)值。

步5 讀取步3和步4同一時(shí)刻2個(gè)波場(chǎng)值，利用互相關(guān)成像條件進(jìn)行成像，從Tmax到0時(shí)刻。

步6 循環(huán)計(jì)算多炮數(shù)據(jù)到Nshot，完成多炮成像結(jié)果疊加。

步7 采用濾波法對(duì)成像結(jié)果進(jìn)行壓制低頻噪聲。

步8 輸出成像結(jié)果。

1.2 CPU/GPU協(xié)同加速逆時(shí)偏移流程

將CUDA編程思想引入到逆時(shí)偏移波場(chǎng)延拓和相關(guān)成像算法中，采用基于CUDA架構(gòu)的CPU/GPU協(xié)同工作模式來(lái)處理逆時(shí)偏移計(jì)算任務(wù)。其主要思想是GPU端負(fù)責(zé)完成波場(chǎng)延拓和相關(guān)成像等并行計(jì)算任務(wù)，CPU端負(fù)責(zé)數(shù)據(jù)輸入和輸出等相關(guān)操作控制工作。

圖2 CPU/GPU協(xié)同加速逆時(shí)偏移流程

首先，將激發(fā)點(diǎn)初始波場(chǎng)值和檢波點(diǎn)最大時(shí)刻波場(chǎng)值由內(nèi)部存儲(chǔ)器(內(nèi)存)傳至設(shè)備存儲(chǔ)器(顯存)中[8]。其次，設(shè)置GPU計(jì)算維度，將GPU多核處理器劃分為相應(yīng)個(gè)數(shù)的計(jì)算塊Block；同時(shí)將每個(gè)Block劃分為若干個(gè)計(jì)算線程Thread；開(kāi)始執(zhí)行炮點(diǎn)波場(chǎng)正傳和檢波點(diǎn)波場(chǎng)反傳計(jì)算任務(wù)。根據(jù)波場(chǎng)計(jì)算相關(guān)成像條件，需要同一時(shí)刻正向傳播的震源波場(chǎng)和反向傳播的檢波點(diǎn)波場(chǎng)，因此互相關(guān)成像計(jì)算也是由GPU來(lái)完成[9，10]，具體實(shí)現(xiàn)流程如圖2所示。

2 系統(tǒng)流程與功能設(shè)計(jì)

2.1 系統(tǒng)設(shè)計(jì)流程

根據(jù)逆時(shí)偏移處理流程，將逆時(shí)偏移系統(tǒng)數(shù)據(jù)處理分為工區(qū)創(chuàng)建、數(shù)據(jù)加載、參數(shù)設(shè)置、作業(yè)發(fā)送和任務(wù)執(zhí)行等多個(gè)步驟，并結(jié)合CPU/GPU協(xié)同并行計(jì)算模式[11]設(shè)計(jì)了逆時(shí)偏移數(shù)據(jù)處理系統(tǒng)流程(見(jiàn)圖3)。

圖3 逆時(shí)偏移系統(tǒng)處理流程

2.2 系統(tǒng)功能設(shè)計(jì)

根據(jù)逆時(shí)偏移的實(shí)際業(yè)務(wù)需求，設(shè)計(jì)并開(kāi)發(fā)了疊前逆時(shí)偏移成像系統(tǒng)，主要完成二維和三維疊前逆時(shí)深度偏移成像功能。系統(tǒng)主要包括資源管理、參數(shù)設(shè)置、數(shù)據(jù)處理、任務(wù)管理、偏移成像和用戶管理等6個(gè)功能模塊(見(jiàn)圖4)。

1)資源管理。對(duì)所有資源節(jié)點(diǎn)進(jìn)行實(shí)時(shí)有效地管理，能夠?qū)崟r(shí)監(jiān)控各個(gè)計(jì)算節(jié)點(diǎn)的運(yùn)行狀態(tài)，獲取節(jié)點(diǎn)的性能參數(shù)，完成計(jì)算節(jié)點(diǎn)的獲取、添加和刪除等操作。

2)參數(shù)設(shè)置。完成逆時(shí)偏移成像前所需要的各類參數(shù)設(shè)置，包括工區(qū)創(chuàng)建、偏移成像前基本參數(shù)設(shè)置和高級(jí)參數(shù)設(shè)置。

圖4 疊前逆時(shí)偏移系統(tǒng)功能模塊圖

3)數(shù)據(jù)處理。完成偏移成像前數(shù)據(jù)的預(yù)處理，包括導(dǎo)入原始數(shù)據(jù)和導(dǎo)入速度模型等操作。

4)偏移成像。包括地震波場(chǎng)正演計(jì)算、檢波點(diǎn)波場(chǎng)反傳計(jì)算和互相關(guān)成像等操作。包括二維數(shù)據(jù)疊前逆時(shí)深度偏移成像和三維數(shù)據(jù)疊前逆時(shí)深度偏移成像功能，其中包括2D CPU RTM、2D GPU RTM、3D CPU RTM和3D GPU RTM。

5)任務(wù)管理。對(duì)逆時(shí)偏移成像任務(wù)的執(zhí)行情況進(jìn)行實(shí)時(shí)監(jiān)測(cè)和管理，能夠監(jiān)管整體任務(wù)和各個(gè)計(jì)算節(jié)點(diǎn)任務(wù)的執(zhí)行情況。包括任務(wù)添加、任務(wù)刪除、任務(wù)重啟、任務(wù)中斷以及查看任務(wù)執(zhí)行時(shí)間等操作，同時(shí)，也能實(shí)時(shí)監(jiān)測(cè)出每個(gè)節(jié)點(diǎn)運(yùn)行性能參數(shù)。

6)用戶管理。主要包括用戶管理、權(quán)限設(shè)置和日志查詢等功能。用戶管理包括用戶的添加、刪除和修改等操作，能夠根據(jù)用戶使用系統(tǒng)功能分配相應(yīng)權(quán)限，根據(jù)權(quán)限不同，用戶登錄后顯示的操作界面也不同，并能夠記錄用戶操作日志。

3 系統(tǒng)主要功能實(shí)現(xiàn)

根據(jù)疊前逆時(shí)偏移數(shù)據(jù)處理流程與實(shí)現(xiàn)算法，通過(guò)研究可視化平臺(tái)開(kāi)發(fā)技術(shù)和疊前逆時(shí)偏移方法，在Linux環(huán)境下應(yīng)用QT圖形界面開(kāi)發(fā)平臺(tái)[12]，設(shè)計(jì)并實(shí)現(xiàn)了一套疊前逆時(shí)偏移成像系統(tǒng)。針對(duì)逆時(shí)偏移實(shí)現(xiàn)過(guò)程中二維CPU、二維GPU、三維CPU和三維GPU偏移算法和處理方法不相同，分別設(shè)計(jì)了二維 CPU逆時(shí)偏移、二維GPU逆時(shí)偏移、三維 CPU逆時(shí)偏移和三維 GPU逆時(shí)偏移等功能。

3.1 疊前逆時(shí)偏移成像實(shí)現(xiàn)

疊前逆時(shí)偏移成像主要包括2D CPU逆時(shí)偏移、2D GPU逆時(shí)偏移、3D CPU逆時(shí)偏移和3D GPU逆時(shí)偏移4個(gè)模塊(見(jiàn)圖5)。筆者以2D GPU逆時(shí)偏移模塊為例描述系統(tǒng)的實(shí)現(xiàn)過(guò)程，其他模塊操作類似：首先創(chuàng)建工區(qū)，導(dǎo)入炮集數(shù)據(jù)和速度模型文件，設(shè)置偏移成像的相關(guān)參數(shù)；選擇2D GPU RTM處理方式和CPU計(jì)算節(jié)點(diǎn)，完成整體偏移成像。2D GPU逆時(shí)偏移具體操作如圖6所示。

圖5 2D CPU逆時(shí)偏移模塊

圖6 2D CPU逆時(shí)偏移成像

偏移成像界面中可以通過(guò)Edit按鈕修改參數(shù)，Inport按鈕載入?yún)?shù)卡文件，參數(shù)設(shè)定后，通過(guò)Export輸出參數(shù)卡文件。所有參數(shù)設(shè)置完成后，點(diǎn)擊Start按鈕發(fā)送作業(yè)開(kāi)始執(zhí)行逆時(shí)偏移計(jì)算，作業(yè)執(zhí)行過(guò)程中點(diǎn)擊Refresh按鈕可以更新進(jìn)度條，查看作業(yè)執(zhí)行進(jìn)度。

3.2 多節(jié)點(diǎn)任務(wù)管理的實(shí)現(xiàn)

圖7 多節(jié)點(diǎn)作業(yè)流程圖

系統(tǒng)設(shè)計(jì)了多節(jié)點(diǎn)間任務(wù)監(jiān)控功能，支持異構(gòu)節(jié)點(diǎn)間作業(yè)的發(fā)送與終止。通過(guò)對(duì)節(jié)點(diǎn)文件信息的讀取，識(shí)別節(jié)點(diǎn)設(shè)備名稱；根據(jù)節(jié)點(diǎn)名與登錄用戶名信息生成唯一的shell文件，作業(yè)發(fā)送或終止命令發(fā)出時(shí)，各節(jié)點(diǎn)執(zhí)行自身匹配的shell文件完成作業(yè)操作。多節(jié)點(diǎn)作業(yè)發(fā)送流程如圖7所示。

在完成多節(jié)點(diǎn)作業(yè)發(fā)送功能基礎(chǔ)上，添加了多節(jié)點(diǎn)作業(yè)進(jìn)度監(jiān)控功能。在多節(jié)點(diǎn)任務(wù)執(zhí)行過(guò)程中，作業(yè)流會(huì)定期更新進(jìn)度文件并保存在工區(qū)指定目錄中。通過(guò)定期發(fā)送更新指令，系統(tǒng)從進(jìn)度文件中讀取作業(yè)進(jìn)度信息，完成進(jìn)度顯示模塊的更新。

4 系統(tǒng)應(yīng)用與效果分析

4.1 系統(tǒng)應(yīng)用可行性分析

為驗(yàn)證系統(tǒng)的可行性及算法集成后的有效性，試驗(yàn)中選取一個(gè)斷層模型的200炮數(shù)據(jù)完成正演計(jì)算和偏移成像試驗(yàn)。試驗(yàn)節(jié)點(diǎn)為1個(gè)Intel E5-2630六核CPU處理器和Nvidia Tesla K10 GPU卡，試驗(yàn)具體參數(shù)如表1所示。

表1 系統(tǒng)軟件測(cè)試

逆時(shí)偏移過(guò)程中采用時(shí)間二階、空間十階的有限差分?jǐn)?shù)值算法完成地震波場(chǎng)正傳和檢波點(diǎn)波場(chǎng)反傳計(jì)算，采用震源歸一化互相關(guān)成像條件提取成像點(diǎn)，并采用拉普拉斯算子濾波法對(duì)圖像進(jìn)行去噪處理，得到最終成像效果。試驗(yàn)中的速度模型、逆時(shí)偏移成像后剖面和采用拉普拉斯算法去噪后剖面分別見(jiàn)圖8(a)、(b)和(c)所示。

圖8 逆時(shí)偏移系統(tǒng)測(cè)試

通過(guò)試驗(yàn)成像結(jié)果發(fā)現(xiàn)，系統(tǒng)成像結(jié)果與源代碼方式成像結(jié)果基本一樣，沒(méi)有明顯區(qū)別。同時(shí)，試驗(yàn)中源代碼計(jì)算與系統(tǒng)計(jì)算的執(zhí)行時(shí)間分別為134.28s和133.46s，系統(tǒng)執(zhí)行時(shí)間略少?？梢?jiàn)，系統(tǒng)完成逆時(shí)偏移從成像效果和執(zhí)行時(shí)間等方面具有可行性，成像效果良好。

4.2 系統(tǒng)2D RTM成像效果分析

表2 試驗(yàn)參數(shù)

為驗(yàn)證系統(tǒng)的成像效果，選取二維鹽丘(Salt-2D)模型進(jìn)行疊前逆時(shí)偏移試驗(yàn)。試驗(yàn)環(huán)境為4個(gè)計(jì)算節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)配備1個(gè)Intel E5-2630六核CPU處理器和Nvidia Tesla K10 GPU卡，試驗(yàn)參數(shù)如表2所示。

在Salt-2D模型試驗(yàn)中，采用時(shí)間二階、空間十階的有限差分?jǐn)?shù)值算法完成波場(chǎng)數(shù)值計(jì)算，采用互相關(guān)成像條件提取成像點(diǎn)，成像后采用拉普拉斯算子濾波法去噪。速度模型如圖9所示，逆時(shí)偏移成像結(jié)果如圖10所示。

從試驗(yàn)結(jié)果可以看出，筆者設(shè)計(jì)的逆時(shí)偏移系統(tǒng)能夠有效的對(duì)地下結(jié)構(gòu)精確成像，深層構(gòu)造成像清楚，整體趨勢(shì)與速度模型數(shù)據(jù)基本一致，地層分界明顯，微幅構(gòu)造成像清晰。

5 結(jié)語(yǔ)

筆者主要圍繞逆時(shí)偏移系統(tǒng)設(shè)計(jì)與應(yīng)用展開(kāi)研究，通過(guò)研究逆時(shí)偏移成像算法和CPU/GPU并行加速技術(shù)，設(shè)計(jì)逆時(shí)偏移系統(tǒng)數(shù)據(jù)處理流程，并根據(jù)CPU和GPU程序代碼、處理方法不同設(shè)計(jì)了2D CPU逆時(shí)偏移、2D GPU逆時(shí)偏移、3D CPU逆時(shí)偏移和3D GPU逆時(shí)偏移等功能模塊，利用QT環(huán)境開(kāi)發(fā)了基于CPU/GPU異構(gòu)平臺(tái)的疊前逆時(shí)偏移成像系統(tǒng)。試驗(yàn)結(jié)果表明系統(tǒng)具有可行性，與源代碼方式相比成像結(jié)果一致，能夠?qū)?fù)雜地下構(gòu)造獲得良好的成像效果。同時(shí)，利用GPU多線程實(shí)現(xiàn)了大規(guī)模并行計(jì)算，大大提高了偏移計(jì)算效率。

圖9 Salt-2D速度模型

圖10 Salt-2D逆時(shí)偏移結(jié)果

[1]陳可洋. 地震波逆時(shí)偏移方法研究綜述[J]. 勘探地球物理進(jìn)展，2010，33(3)：153～159.

[2] Liu Y，Chang X，Jin D，et al. Reverse time migration of multiples for subsalt imaging[J]. Geophysics，2011，76(5)：209～216.

[3] Costa J C，Silva Neto F A.Obliquity-correction imaging condition for reverse time migration[J]. Geophysics，2009，74(3)：57～66.

[4] 高新成，李春生.基于CPU/GPU協(xié)同加速疊前逆時(shí)偏移方法研究[J].陜西理工學(xué)院學(xué)報(bào)(自然科學(xué)版)，2016(2)：44～49.

[5] 石穎，柯璇，張瑩瑩. 逆時(shí)偏移邊界條件與存儲(chǔ)策略分析[J].地球物理學(xué)進(jìn)展，2015(2)：581～585.

[6] Clapp R G. Reverse time migration with random boundaries[J]. 79th Annual International Meeting，SEG Expanded Abstracts，2009，28：2809～2813.

[7]Jiaze He，F(xiàn)uh-Gwo Yuan.Damage identification for composite structures using a cross-correlation reverse-time migration technique[J]. Structural Health Monitoring，2015， 146：224～226.

[8] Ying Shi，Yanghua Wang. Reverse time migration of 3D vertical seismic profile data[J]. Geophysics，2016，81(1)：31～38.

[9] 石穎，陸加敏，柯璇，等. 基于GPU并行加速的疊前逆時(shí)偏移方法[J].東北石油大學(xué)學(xué)報(bào)，2012，36(4)：111～115.

[10] 劉守偉，王華忠，陳生昌，等. 三維逆時(shí)偏移GPU/CPU機(jī)群實(shí)現(xiàn)方案研究[J]. 地球物理學(xué)報(bào)，2013，56(10)：3487～3496.

[11] 王維紅，柯璇，郭雪豹，等. 基于QT的Walkaway VSP逆時(shí)偏移軟件開(kāi)發(fā)研究[J]. 石油物探，2015，54(4)：452～458.

[12] 陳楠，祝媛媛，張光德，等. 基于QT的地震勘探可擴(kuò)展平臺(tái)研發(fā)與應(yīng)用[J]. 華北地震科學(xué)，2013，32(4)：31～36.

[編輯] 辛長(zhǎng)靜

2016-09-18

國(guó)家自然科學(xué)基金項(xiàng)目(41574117)；黑龍江省教育科學(xué)規(guī)劃重點(diǎn)課題(GJB1215013)

高新成(1979-)，男，博士，副教授，現(xiàn)主要從事地震數(shù)據(jù)處理與智能計(jì)算方面的教學(xué)與研究工作；E-mail：gxc@nepu.edu.cn。

TP399

A

1673-1409(2016)34-0015-06

[引著格式]高新成，石穎.基于CPU/GPU異構(gòu)平臺(tái)的疊前逆時(shí)偏移成像系統(tǒng)[J].長(zhǎng)江大學(xué)學(xué)報(bào)(自科版)，2016,13(34)：15～20.