基于gpu加速和隨機(jī)邊界的三維彈性波逆時偏移算法

基于gpu加速和隨機(jī)邊界的三維彈性波逆時偏移算法

1.gpu并行計算技術(shù)逆時偏移需要精確地求解雙向聲波或彈性波的波動方程，反向地向外推地震數(shù)據(jù)，并選擇合適的成像條件來精確成像。目前逆時偏移是公認(rèn)的最能精確成像的偏移技術(shù),原因是其適用于任意復(fù)雜介質(zhì)、不存在任何傾角限制、能夠?qū)θ我忸愋偷牟ㄟM(jìn)行數(shù)值模擬成像相比CPU,GPU旨在增強計算機(jī)圖形處理顯示能力,天生具備高并行性和高密集度數(shù)值計算的巨大潛能。目前,GPU硬件平臺和CUDA編程日漸成熟,越來越多的業(yè)內(nèi)人士開始尋求利用GPU超強的并行計算能力解決三維逆時偏移、三維模型數(shù)值模擬、速度建模以及全波形反演等地球物理勘探的前沿問題。目前將GPU并行計算應(yīng)用于海量地球物理數(shù)據(jù)運算已取得大規(guī)模提速,為大規(guī)模實際地球物理資料的逆時偏移實現(xiàn)搭建了平臺基礎(chǔ)。針對逆時偏移需要巨大臨時文件存儲量的問題,業(yè)內(nèi)提出了Checkpointing方法、邊界存儲策略等方法來降低硬盤存儲,但仍需存儲相當(dāng)大的數(shù)據(jù)。貫徹“用計算換存儲”策略的隨機(jī)邊界存儲方法目前三維逆時偏移仍受限于計算量和存儲量巨大等問題而難以推廣,基于GPU并行加速和隨機(jī)邊界存儲策略的三維彈性波逆時偏移算法研究,提高逆時偏移效率,是一個有很強現(xiàn)實意義的課題。2.三維彈性波方程的逆偏移原理2.1d交錯網(wǎng)格數(shù)值離散3D各向同性介質(zhì)中的一階速度—應(yīng)力彈性波方程如下:式(1)中,x、y、z表示空間直角空間坐標(biāo)方向,t代表時間,ρ代表密度。v本文采用3D交錯網(wǎng)格對方程(1)進(jìn)行數(shù)值離散,交錯網(wǎng)格中速度分量v表1、圖1展示了3D交錯網(wǎng)格中各彈性波波場分量、彈性參數(shù)的空間分布。在3D交錯網(wǎng)格空間中對式(1)進(jìn)行有限差分,得到時間2階,空間上任意偶數(shù)階的有限差分格式,以式(1)中第1個方程舉例:式(2)中,Δx、Δy、Δz代表空間x、y、z方向離散間隔;Δt代表時間離散間隔;i、j、k、n表示空間x、y、z方向和時間上的離散點序號；i2.2受內(nèi)部彈性介質(zhì)的傳播方向本文通過傅立葉變換和矩陣特征值分析可得穩(wěn)定性條件在三維各向同性介質(zhì)中,彈性波的傳播速度是與傳播方向無關(guān)的,因此x,y,z三方向應(yīng)滿足穩(wěn)定條件是相同的。2.3隨機(jī)速度隨機(jī)擴(kuò)散波場模擬本文利用隨機(jī)邊界存儲策略進(jìn)行逆時偏移,基本思路是:選擇合適的隨機(jī)方程,生成的隨機(jī)邊界,消除彈性波方程延拓中人工邊界和自由邊界條件下反射波的相干性,使得邊界反射不能進(jìn)行有效成像。具體實現(xiàn)實現(xiàn)過程如下:將有限計算空間向外擴(kuò)充一定的區(qū)域,在外擴(kuò)的區(qū)域內(nèi)填充隨機(jī)速度,當(dāng)彈性波場越過模型邊界向外繼續(xù)傳播到隨機(jī)速度區(qū)域時,內(nèi)充隨機(jī)速度對波前面進(jìn)行隨機(jī)化,使波場信息變成隨機(jī)噪聲,反傳回真實的速度區(qū)域,破壞模型邊界反射的相干性,進(jìn)而使模型邊界反射不能進(jìn)行有效成像。注意:隨機(jī)邊界機(jī)速度也要滿足穩(wěn)定性條件公式(3)。本文使用3D各向同性的地質(zhì)模型來驗證隨機(jī)邊界的效果。三維均勻介質(zhì)模型大小為500m×500m×500m,采樣網(wǎng)格Δx、Δy、Δz為5m,縱波速度v2.4矢量橫波標(biāo)量化據(jù)馬德堂等,縱波為無旋場,橫波為無散場該方法所得縱波為標(biāo)量,橫波則為矢量,需要進(jìn)行矢量橫波標(biāo)量化,這樣獲得的標(biāo)量橫波可以直接用于逆時偏移的成像。2.5成像剖面的選取彈性波逆時偏移的另一個關(guān)鍵是偏移成像條件的選取,這將直接影響成像剖面的最終質(zhì)量。本文采用歸一化互相關(guān)成像條件,通過計算兩個波場的零延遲互相關(guān),并做歸一化處理,進(jìn)行偏移成像。3.gpu并行性分析大型地震資料的逆時偏移因其巨大的計算量和數(shù)據(jù)存儲需求而備受限制,上文分析可得逆時偏移的延拓和成像算法是可分解和可線性疊加的,因而具備大規(guī)模并行計算的可能。待移植程序的計算部分是否具有并行性以及計算量大小是程序移植之前必須考慮的2個問題。若不能同時滿足這兩點要求,即是程序成功移植到GPU平臺,也會出現(xiàn)并行性不高,效率沒有明顯提升的現(xiàn)象,有時甚至?xí)档?。串行計算中三維彈性波逆時延拓由CPU單一線程順序?qū)崿F(xiàn),GPU平臺可通過多線程并行執(zhí)行。GPU內(nèi)沒有(或有少量)高速緩存,事實證明訪問共享存儲器的速度和訪問寄存器的速度相差不多,因此可以通過訪問共享存儲器來實現(xiàn)不同線程間的最小通信延遲。龍桂華等分析發(fā)現(xiàn)GTX295顯卡上利用片內(nèi)共享存儲計算1024×1024兩矩陣乘法時速度大約為利用片內(nèi)共享內(nèi)存時計算速度的4.5倍4.空間大時延拓的隨機(jī)邊界存儲策略逆時偏移所采用的成像條件需要對同一時刻的炮點和檢波點波場進(jìn)行相關(guān)計算,而炮點和檢波點波場是相互反向的過程,不能同時由計算,因此常規(guī)逆時偏移通常選擇在炮點波場的正向延拓時記錄整個波場的值,在檢波點波場的逆時延拓過程中直接讀入所需同一時刻的震源波場值。這種需要存儲整個模型波場值的方法硬盤存儲量是巨大的,大約x×y×z×t。而巨大硬盤存儲要求是制約3D資料彈性波逆時偏移推廣的又一瓶頸。另外,震源波場正向延拓記錄數(shù)值、檢波點波場逆時延拓讀取數(shù)值等操作都將導(dǎo)致I/O訪問過大,耗費大量時間,嚴(yán)重影響計算效率。對于CUDA這種異構(gòu)編程模式而言,CPU/GPU通信是必須考慮的問題,需要最大限度地減少I/O訪問、讀取操作。因此,本文采用隨機(jī)邊界存儲策略,在炮點波場正向延拓過程中,使用隨機(jī)速度構(gòu)建隨機(jī)邊界,只需保存最后t時刻的炮點波場值(包括隨機(jī)邊界范圍在內(nèi)的全部炮點波場值),在檢波點波場進(jìn)行逆時延拓的同時,通過讀取保存的最后t時刻隨機(jī)波場快照和隨機(jī)邊界同步逆時重構(gòu)炮點波場。這樣,相同時刻的炮點波場值和檢波點波場值能隨時計算出來,進(jìn)而進(jìn)行偏移成像。雖然這種方法對炮點波場延拓進(jìn)行了2次計算,但最大限度地降低了硬盤空間存儲量和文件讀寫、拷貝的工作量,相比于GPU極高的并行計算能力,增加的計算量可忽略不計。5.地表剖面位置及三維傾斜模型三維傾斜模型,底部有一高速圓形隆起,模型1500m×1500m×1500m,采樣網(wǎng)格Δx、Δy、Δz為5m,時間間隔為0.5ms,記錄1.5s。一共布設(shè)10條測線,測間距為150m,測線內(nèi)炮間距為100m,每條測線布設(shè)15炮,每炮30×30道接收,處理面元50m×50m,在地表處(z=0)接收。三維傾斜模型的縱橫波速度模型以及逆時偏移剖面如圖3所示。結(jié)果顯示3D彈性波逆時偏移算法能夠比較準(zhǔn)確對模型進(jìn)行成像,識別出其中的水平地層、傾斜地層和底部穹形隆起。本次3D傾斜模型一共正演150炮,逆時偏移150炮。比較一炮數(shù)據(jù)計算時間,CPU算法:6301.62s,GPU并行+隨機(jī)邊界算法(GTX560Ti,核心448):137.41s,加速比達(dá)45.86,計算效率