一種三維各向異性彈性波數(shù)值模擬方法及系統(tǒng)的制作方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及一種三維各向異性彈性波數(shù)值模擬方法及系統(tǒng)，屬于地球物理勘探領(lǐng) 域。
【背景技術(shù)】
[0002] 地震波數(shù)值模擬以地震波在地下介質(zhì)中的傳播理論為基礎(chǔ)，在勘探地震學(xué)與天然 地震學(xué)中得到了廣泛的應(yīng)用。目前常規(guī)的三維聲波方程、彈性波各向同性數(shù)值模擬已經(jīng)廣 泛的應(yīng)用于數(shù)值模擬、成像、以及反演等各個(gè)地球物理領(lǐng)域中。但是對(duì)于如何有效的實(shí)現(xiàn)大 尺度采集方式(例如寬方位采集)和復(fù)雜各向異性介質(zhì)（例如水平橫向各向同性HTI或者是 正交各向異性）的三維彈性波波場(chǎng)數(shù)值模擬研究仍然存在很大的挑戰(zhàn)，在實(shí)際應(yīng)用中未得 到廣泛的使用。此外常規(guī)的基于CPU的大規(guī)模三維聲波、彈性波數(shù)值模擬，通常需要大量專 用集群計(jì)算資源。無(wú)論從硬件成本還是計(jì)算能耗上來(lái)說(shuō)其成本都十分昂貴。因此快速的提 高計(jì)算性能、顯著的降低計(jì)算成本對(duì)于實(shí)現(xiàn)三維彈性波各向異性數(shù)值模擬的實(shí)際應(yīng)用具有 重要意義。
[0003] 近十年中，利用圖形處理單元(GPU)進(jìn)行計(jì)算密集型應(yīng)用的加速實(shí)現(xiàn)已經(jīng)得到了 突飛猛進(jìn)式的發(fā)展。圖形處理單元(GPU)由于其具有高速的內(nèi)存帶寬，相較于CPU至少高出 兩個(gè)數(shù)量級(jí)的計(jì)算處理核心，更適合并行計(jì)算的單指令多數(shù)據(jù)(SMD)計(jì)算模式，以及更低 的能耗成本，正廣泛的應(yīng)用到計(jì)算科學(xué)的相關(guān)領(lǐng)域。對(duì)于勘探地球物理領(lǐng)域，對(duì)于使用圖形 處理單元(GPU)的興趣也在顯著增強(qiáng)，越來(lái)越多的研究已經(jīng)將GPU用于加速地震處理中的核 心算法，例如地震數(shù)值模擬、地震成像、地震高精度反演等。

【發(fā)明內(nèi)容】

[0004] 經(jīng)過(guò)研究發(fā)現(xiàn)，現(xiàn)有技術(shù)及其存在以下問(wèn)題：
[0005] 將時(shí)間域有限差分方法應(yīng)用到GHJ設(shè)備上實(shí)現(xiàn)復(fù)雜介質(zhì)波傳播算法上的研究很 少，尤其對(duì)于如何處理存在大量多節(jié)點(diǎn)數(shù)據(jù)交換的大規(guī)模三維問(wèn)題研究就更少(Heinrich et al，2014)〇
[0006] 目前存在一些GPU地震波傳播模擬的實(shí)現(xiàn)方案，對(duì)于三維GPU集群的實(shí)現(xiàn)方案我們 以（龍桂華等）的方案為例進(jìn)行描述。該方案利用區(qū)域分解技術(shù)將單個(gè)GPU上不能計(jì)算的地 質(zhì)體沿Z軸方向進(jìn)行粗粒度分解，采用消息傳遞接口（MPI)交換邊界數(shù)據(jù)，從而運(yùn)用MPI + CUDA的方式實(shí)現(xiàn)了大尺度三維地震波場(chǎng)的數(shù)值模擬。但是該方法存在很大的一個(gè)問(wèn)題是， 采用GHJ集群計(jì)算的加速比相較于單個(gè)GPU與CPU來(lái)說(shuō)顯著下降，造成這一結(jié)果的原因是GPU 計(jì)算的大部分耗時(shí)消耗在GPU集群節(jié)點(diǎn)間的GPU到CPU以及CPU到GPU設(shè)備之間的數(shù)據(jù)拷貝 上。
[0007] 本發(fā)明所要解決的技術(shù)問(wèn)題是，針對(duì)現(xiàn)有技術(shù)的不足，提供一種將GPUDirect技術(shù) 與有限差分?jǐn)?shù)值模擬方法結(jié)合得到的基于GPU Direct優(yōu)化通信的三維各向異性彈性波數(shù) 值模擬方法及系統(tǒng)，為大規(guī)模數(shù)值模擬提供有力保證。
[0008] 本發(fā)明解決上述技術(shù)問(wèn)題的技術(shù)方案如下:一種三維各向異性彈性波數(shù)值模擬方 法，具體包括以下步驟：
[0009] 步驟1:建立介質(zhì)模型，對(duì)介質(zhì)模型進(jìn)行網(wǎng)格離散得到多個(gè)網(wǎng)格點(diǎn)；
[0010] 步驟2:計(jì)算震源函數(shù)，根據(jù)震源函數(shù)計(jì)算每個(gè)網(wǎng)格點(diǎn)上的壓力值；
[0011] 步驟3:將三維各向異性彈性波方程轉(zhuǎn)換為傳播方程，將每個(gè)網(wǎng)格點(diǎn)上的壓力值帶 入傳播方程進(jìn)行計(jì)算，得到每一刻的波場(chǎng)值；
[0012] 步驟4:根據(jù)波場(chǎng)值確定每個(gè)網(wǎng)格點(diǎn)的計(jì)算區(qū)域，進(jìn)行分區(qū)并對(duì)分區(qū)邊界數(shù)據(jù)進(jìn)行 數(shù)據(jù)交換，完成彈性波數(shù)值模擬。
[0013] 本發(fā)明的有益效果是:從應(yīng)力應(yīng)變方程出發(fā)，實(shí)現(xiàn)了利用圖形處理單元(GPU)加速 復(fù)雜介質(zhì)下的三維彈性波數(shù)值模擬，并針對(duì)三維問(wèn)題引入GPU設(shè)備所面對(duì)的由于區(qū)域分解 所造成的多節(jié)點(diǎn)間、節(jié)點(diǎn)內(nèi)通信瓶頸問(wèn)題進(jìn)行了深入的研究和分析，提出了利用GPU Direct技術(shù)加速數(shù)據(jù)傳輸?shù)膶?shí)現(xiàn)方案，避免了大量的CPU到GPU，GPU到CPU的數(shù)據(jù)拷貝，實(shí)現(xiàn) 了優(yōu)化通信瓶頸的問(wèn)題。通過(guò)GPU圖形加速設(shè)備并引入GPU Direct通信優(yōu)化策略，可以顯著 的提升整體計(jì)算性能，可以用更低的硬件成本和更少的時(shí)間實(shí)現(xiàn)的三維各向異性彈性數(shù)值 模擬，為各種依賴于波動(dòng)方程正演模擬的算法如逆時(shí)偏移，全波形反演的應(yīng)用提供有力的 保證。
[0014] 在上述技術(shù)方案的基礎(chǔ)上，本發(fā)明還可以做如下改進(jìn)。
[0015] 進(jìn)一步，所述步驟4具體包括以下步驟：
[0016] 步驟4.1:根據(jù)所有波場(chǎng)值確定每個(gè)網(wǎng)格點(diǎn)的計(jì)算區(qū)域，并采用吸收邊界的方式確 定計(jì)算區(qū)域邊界，在計(jì)算區(qū)域內(nèi)模擬地下介質(zhì)中波的傳播；
[0017] 步驟4.2:對(duì)所有計(jì)算區(qū)域進(jìn)行分區(qū)，對(duì)相鄰分區(qū)的邊界數(shù)據(jù)進(jìn)行數(shù)據(jù)交換，得到 模擬彈性波數(shù)據(jù)，完成彈性波數(shù)值模擬。
[0018] 進(jìn)一步，所述步驟4.2中邊界數(shù)據(jù)的采用GPU-Direct技術(shù)進(jìn)行數(shù)據(jù)交換。
[0019] 采用上述進(jìn)一步方案的有益效果是，利用GPU加速三維彈性波數(shù)值模擬，顯著提高 計(jì)算速度，相較于常規(guī)CHJ集群實(shí)現(xiàn)，可以達(dá)到20-30倍的速度提升。
[0020] 進(jìn)一步，所述步驟4.2中的數(shù)據(jù)交換具體包括：
[0021] 沿介質(zhì)模型中計(jì)算區(qū)域邊界數(shù)據(jù)變化最慢的軸向進(jìn)行區(qū)域分解，得到多個(gè)分區(qū)， 將所有分區(qū)分配到多個(gè)GPU中，每個(gè)分區(qū)獨(dú)立在一個(gè)GPU上執(zhí)行計(jì)算;每?jī)蓚€(gè)相鄰分區(qū)的邊 界數(shù)據(jù)再進(jìn)行交換。
[0022]采用上述進(jìn)一步方案的有益效果是，利用GPU Direct技術(shù)加速節(jié)點(diǎn)間、節(jié)點(diǎn)內(nèi)數(shù) 據(jù)傳輸?shù)膶?shí)現(xiàn)方案，我們的方案因其避免了大量的CPU到GPU，GPU到CPU的數(shù)據(jù)拷貝，從而實(shí) 現(xiàn)了優(yōu)化通信瓶頸的問(wèn)題。
[0023] 進(jìn)一步，所述步驟1具體包括以下步驟：
[0024] 步驟1.1:根據(jù)地質(zhì)背景條件、實(shí)際測(cè)得的巖石物理測(cè)試數(shù)據(jù)和測(cè)井資料數(shù)據(jù)建立 介質(zhì)模型；
[0025]步驟1.2:采用形狀規(guī)則的三維網(wǎng)格對(duì)介質(zhì)模型進(jìn)行網(wǎng)格離散，得到網(wǎng)格點(diǎn)。
[0026]采用上述進(jìn)一步方案的有益效果是，將介質(zhì)模型離散成網(wǎng)格點(diǎn)，網(wǎng)格點(diǎn)作為波場(chǎng) 值的載體，由三維結(jié)構(gòu)的多個(gè)網(wǎng)格點(diǎn)構(gòu)成波場(chǎng)，作為波場(chǎng)傳播的初始條件。
[0027]進(jìn)一步，所述震源函數(shù)在空間上采用高斯函數(shù)，在時(shí)間上采用Ricker子波，所述震 源函數(shù)的具體公式為：
[0028] s(x，y，z，t)=g(x，y，z) · f(t)公式（1)
[0029] 其中，以〇 =(卜2(31侃)2)6叉？（-(31侃）2)公式（2)
[0030]
[0031 ] 式中：t表示時(shí)間，fo表示Ricker子波的中心頻率，模型計(jì)算中fo = 15Ηζ，β為常數(shù)； (義0,7〇,2())為震源的空間位置，1、7和2分別為1軸、7軸和2軸方向上的位置。
[0032] 進(jìn)一步，所述步驟3具體包括以下步驟：
[0033]步驟3.1:將三維各向異性彈性波方程中的微分用差分近似替代，得到相應(yīng)的有限 差分格式的傳播方程，所述傳播方程中空間采樣步長(zhǎng)和時(shí)間采樣步長(zhǎng)滿足該數(shù)值格式的穩(wěn) 定性條件；
[0034]步驟3.2:采用區(qū)域分解的方式將每個(gè)網(wǎng)格點(diǎn)上的壓力值帶入傳播方程進(jìn)行計(jì)算， 得到每一刻的波場(chǎng)值。
[0035]本發(fā)明解決上述技術(shù)問(wèn)題的技術(shù)方案如下:一種三維各向異性彈性波數(shù)值模擬系 統(tǒng)，包括建模模塊、震源模塊、傳播模塊和數(shù)據(jù)交換模塊；
[0036] 所述建模模塊用于建立介質(zhì)模型，對(duì)介質(zhì)模型進(jìn)行網(wǎng)格離散得到多個(gè)網(wǎng)格點(diǎn)；
[0037] 所述震源模塊用于計(jì)算震源函數(shù)，根據(jù)震源函數(shù)計(jì)算每個(gè)網(wǎng)格點(diǎn)上的壓力值； [0038]所述傳播模塊用于將三維各向異性彈性波方程轉(zhuǎn)換為傳播方程，將每個(gè)網(wǎng)格點(diǎn)上 的壓力值帶入傳播方程進(jìn)行計(jì)算，得到每一刻的波場(chǎng)值；
[0039]所述數(shù)據(jù)交換模塊用于根據(jù)波場(chǎng)值確定每個(gè)網(wǎng)格點(diǎn)的計(jì)算區(qū)域，進(jìn)行分區(qū)并對(duì)分 區(qū)邊界數(shù)