vi介質(zhì)高精度疊前逆時(shí)偏移方法研究

vi介質(zhì)高精度疊前逆時(shí)偏移方法研究

0vti介質(zhì)四階波動(dòng)方程重疊前的逆偏移（rtm）是一種高精度的地震學(xué)技術(shù)，不僅可以打破地層平面的界限，還可以成像各種波，如旋轉(zhuǎn)波、角度波和多個(gè)反射波。20世紀(jì)80年代，許多科學(xué)家提出了這種算法。然而，由于rtm方法對(duì)硬件的要求很高，它一直在發(fā)展緩慢。近年來，隨著計(jì)算機(jī)硬件技術(shù)的發(fā)展，random反向偏移法越來越受到重視。在地球上的各向異性介質(zhì)中廣泛存在。為了在各向異性介質(zhì)中采用各向異性介質(zhì)模型進(jìn)行逆偏移，成像精度通常降低。例如，反射波的返回不允許，能量也不聚焦。為了提高反向偏移技術(shù)的精度和效率，有必要考慮不同向異性介質(zhì)的影響。公共時(shí)間偏移的估計(jì)可以通過波場合成。然而，不僅各向異性介質(zhì)的彈性波方程相對(duì)復(fù)雜，而且很難完全分離垂直波場，計(jì)算成本高，這無疑對(duì)高成本的反向偏移提出了更高的計(jì)算要求。為了避免彈性波方程的計(jì)算，使用每個(gè)向異性介質(zhì)的準(zhǔn)p波方程。垂直向同性介質(zhì)（vti）是每個(gè)向異性介質(zhì)的一般模型。因此，研究第一部分的逆p波方程具有重要的意義，在其他各向異性介質(zhì)的高精度逆偏移中具有重要意義。目前,存儲(chǔ)量大、計(jì)算成本高以及低頻噪音仍是限制逆時(shí)偏移技術(shù)發(fā)展的主要問題.近年來,多名學(xué)者提出或改進(jìn)了RTM的存儲(chǔ)方案,Symes(Symes,2007)提出了用checkpointing技術(shù)來降低波場存儲(chǔ)量的思想;Dussaud(Dussaudetal.,2008)等比較了幾種常用的存儲(chǔ)策略,認(rèn)為checkpointing技術(shù)降低了波場存儲(chǔ)量的同時(shí)較明顯增加了計(jì)算量;Clapp(Clapp,2008,2009)提出了隨機(jī)邊界存儲(chǔ)方法,在很大程度上節(jié)約了存儲(chǔ)量.王保利(王保利等,2012)等詳細(xì)比較了幾種常用的邊界存儲(chǔ)策略,提出了有效邊界存儲(chǔ)策略.胡昊(胡昊等,2013)對(duì)比分析了隨機(jī)邊界和單層吸收邊界,認(rèn)為隨機(jī)邊界會(huì)出現(xiàn)邊界漫反射,單層吸收邊界會(huì)損失部分能量等.康瑋(康瑋和程玖兵,2012)分析了疊前逆時(shí)偏移低頻噪音干擾去除的一些常用方法,認(rèn)為基于GPU平臺(tái)的逆時(shí)偏移,采用Laplace去噪是目前最經(jīng)濟(jì)的選擇.最近幾年圖形處理器(GraphicProcessingUnit)的迅猛發(fā)展為浮點(diǎn)數(shù)據(jù)計(jì)算帶來了大幅度的效率提升.劉紅偉(劉紅偉等,2010)和李博(李博等,2010)等借助GPU加速技術(shù)實(shí)現(xiàn)了有限差分逆時(shí)偏移算法,提極大地提高了計(jì)算效率.為解決各向異性介質(zhì)彈性波數(shù)值模擬的問題,Alkhalifah(Alkhalifah,1998,2000)提出了聲波近似思想,給出了VTI介質(zhì)四階波動(dòng)方程和二階方程組.研究表明,聲學(xué)近似對(duì)縱波的運(yùn)動(dòng)學(xué)特征幾乎無影響,同時(shí)避免了彈性波數(shù)值模擬的缺陷.在Alkhalifah研究的基礎(chǔ)上,一些研究人員給出了其他形式的二階偏微分方程組(KliéandToro,2001;Zhouetal.,2006;Hestholm,2007;Duetal.,2008;Duvenecketal.,2008).Hestholm(Hestholm,2009)對(duì)二階VTI方程組進(jìn)行了數(shù)學(xué)變換,組合成為六個(gè)一階耦合偏微分方程.Zhou(Zhouetal.,2006),Fletcher(Fletcheretal.,2009),Fowler(Fowleretal.,2010)等利用聲波近似原理推導(dǎo)出了TTI介質(zhì)擬聲波方程.Bube(Bubeetal.,2012),Zhang(ZhangandZhang,2009;ZhangandSun,2009),Bakker(BakkerandDuveneck,2011)等針對(duì)TTI介質(zhì)波場模擬不穩(wěn)定現(xiàn)象,提出了一些改進(jìn)的思想以及技術(shù)方法.張巖(張巖和吳國忱,2013)回顧了各向異性介質(zhì)逆時(shí)偏移發(fā)展歷史,對(duì)TTI介質(zhì)疊前逆時(shí)偏移成像的研究現(xiàn)狀、進(jìn)展進(jìn)行了概述.沈銘成(沈銘成等,2014)針對(duì)VTI介質(zhì)聲學(xué)近似帶來的偽橫波干擾以及算法穩(wěn)定性問題,提出了VTI介質(zhì)的解耦合聲波近似方程,該方程可以使qp波波場傳播更穩(wěn)定,同時(shí)消除了偽sv波的干擾.程玖兵(程玖兵等,2013,2014)研究了各向異性介質(zhì)qp波傳播模式:偽純模式波動(dòng)方程和分離純模式標(biāo)量波,為各向異性介質(zhì)的分離模式的波場的傳播過程提供了描述工具.交錯(cuò)網(wǎng)格(Virieux,1984)相對(duì)于常規(guī)網(wǎng)格精度高,本文在Hestholm(Hestholm,2009)研究的基礎(chǔ)上,推導(dǎo)出了該方程組炮點(diǎn)和檢波點(diǎn)逆時(shí)延拓的交錯(cuò)網(wǎng)格高階有限差分格式,針對(duì)算法計(jì)算量大、存儲(chǔ)量大的問題,給出了基于GPU平臺(tái)的改進(jìn)有效完全匹配層(PML)吸收邊界條件,不僅提高了計(jì)算效率同時(shí)降低了存儲(chǔ)量.應(yīng)用震源歸一化互相關(guān)成像條件,利用拉普拉斯算子去除低頻噪音,利用Hess2DVTI模型證明本文方法的正確性及有效性.1在vsti介質(zhì)中，第一學(xué)期的準(zhǔn)p波方程的傳輸點(diǎn)和檢測點(diǎn)的逆時(shí)延1.1hestholm表示的一階準(zhǔn)p波方程為二維vsi介質(zhì)的中一階準(zhǔn)p波方程式中:P為應(yīng)力,ρ為密度,V1.2時(shí)間離散擬合利用交錯(cuò)網(wǎng)格(董良國等,2000)對(duì)式(1)進(jìn)行有限差分離散,推導(dǎo)出準(zhǔn)P波的逆時(shí)延拓差分格式(本文只離散式(1-a)、(1-f),同理可對(duì)其他式離散):式中:Δx、Δz、Δt分別為空間x方向離散步長、空間z方向離散步長、時(shí)間離散步長;i、j為空間離散位置坐標(biāo);n為時(shí)間取樣點(diǎn)數(shù);C基于PML邊界條件的分裂思路(王守東,2003),對(duì)式(1)研究區(qū)域四周引入完全匹配層,可以得到完全匹配層控制方程,本文只對(duì)式(1-a)、(1-f)進(jìn)行匹配層引入,同理可得其余式匹配層控制方程為同理對(duì)式(3)進(jìn)行交錯(cuò)網(wǎng)格有限差分格式離散,可以得到完全匹配層正向延拓和逆時(shí)延拓交錯(cuò)網(wǎng)格差分格式.1.3階準(zhǔn)p波交錯(cuò)網(wǎng)格差分格式對(duì)于時(shí)間2階差分精度,空間2N階差分精度的一階準(zhǔn)P波交錯(cuò)網(wǎng)格差分格式,本文采用的穩(wěn)定性條件(董良國等,2000)為式中:Δt,Δx,Δz分別為時(shí)間和空間采樣間隔,max[v2基于gpu加速的pml吸收邊界條件存儲(chǔ)策略逆時(shí)偏移方法包括正傳波場計(jì)算、反傳波場計(jì)算以及二者的零延遲互相關(guān)成像,若保存全部歷史時(shí)刻的正傳波場,會(huì)出現(xiàn)較大的存儲(chǔ)需求.例如:2000×1000網(wǎng)格點(diǎn)數(shù)的速度模型,時(shí)間采樣點(diǎn)數(shù)取10000,存儲(chǔ)數(shù)據(jù)類型是4字節(jié)浮點(diǎn)型,則需要約80Gb來存儲(chǔ)所有時(shí)刻的正傳波場,目前的計(jì)算機(jī)硬件幾乎無法承受.高存儲(chǔ)成本是限制逆時(shí)偏移技術(shù)發(fā)展的主要障礙,計(jì)算機(jī)硬件技術(shù)的發(fā)展在一定程度上可以帶動(dòng)存儲(chǔ)能力的提高,但是,仍需尋求有效降低存儲(chǔ)的逆時(shí)偏移算法.鑒于此,本文提出了改進(jìn)的基于GPU加速的PML吸收邊界存儲(chǔ)策略,其可在降低存儲(chǔ)成本的同時(shí),實(shí)現(xiàn)高效率、高精度的疊前逆時(shí)偏移成像.2.1改進(jìn)的pml邊界存儲(chǔ)策略PML吸收邊界條件使波場在PML吸收層衰減,是目前效果最好的邊界吸收條件.本文針對(duì)交錯(cuò)網(wǎng)格,研究了一種改進(jìn)的PML邊界存儲(chǔ)方法,該方法存儲(chǔ)量小、計(jì)算成本低、吸收效果好.該方法借鑒了隨機(jī)邊界逆向傳播震源波場的思想,保存和利用最后兩個(gè)時(shí)刻的震源波場,有別于隨機(jī)邊界方法的是,文中方法也存儲(chǔ)正向傳播的各時(shí)刻邊界波場.首先,在震源波場正向傳播的過程中把邊界的波場存儲(chǔ)起來,然后,在后兩時(shí)刻的波場逆時(shí)回推前一時(shí)刻的波場的過程中,用存儲(chǔ)的邊界波場代替逆時(shí)回推過程的邊界波場,這樣可以保證在模擬計(jì)算區(qū)域,每個(gè)時(shí)刻波場都具有完整的能量,因此后兩時(shí)刻的波場可以恢復(fù)前一時(shí)刻的波場.王保利(王保利等,2012)提出了有效邊界存儲(chǔ)策略,其方法是修改波場逆向傳播的邊界條件,保存PML區(qū)域內(nèi)的若干層內(nèi)所有時(shí)刻的波場值,利用最后兩個(gè)時(shí)刻波場通過反傳播技術(shù)得到之前任意時(shí)刻的正確的波場.其邊界波場存儲(chǔ)方案如圖1所示:虛線框包括了PML吸收邊界和內(nèi)部模擬計(jì)算區(qū)域,虛線和相鄰的的實(shí)線框所組成的白色區(qū)域?yàn)镻ML吸收邊界層,兩個(gè)實(shí)線框所組成的灰色區(qū)域?yàn)橐涗浀倪吔绮▓?實(shí)線邊框內(nèi)是不包括吸收邊界的計(jì)算區(qū)域.本文提出了一種改進(jìn)方法,如圖2所示,相比于圖1所示的存儲(chǔ)方案,本文改進(jìn)方法所存儲(chǔ)的邊界波場是計(jì)算區(qū)域的波場,而非PML吸收邊界區(qū)域的波場,邊界波場在內(nèi)部模擬計(jì)算區(qū)域,不會(huì)有任何能量衰減,也在減少存儲(chǔ)量的同時(shí),縮短了讀寫邊界波場的時(shí)間.下面將詳細(xì)的介紹改進(jìn)的PML邊界存儲(chǔ)策略以及VTI介質(zhì)一階準(zhǔn)P波的逆時(shí)偏移實(shí)現(xiàn)過程.公式(2)為推導(dǎo)出來的VTI介質(zhì)一階準(zhǔn)P波交錯(cuò)網(wǎng)格差分格式,包括六個(gè)相互耦合的變量:波場P,速度分量v為驗(yàn)證改進(jìn)的PML邊界條件能否正確逆時(shí)回推震源波場,設(shè)計(jì)了200×200網(wǎng)格的均勻各向異性介質(zhì),其中各向異性參數(shù)η為0.5,動(dòng)校正速度和垂向速度均為3000m/s,完全匹配層(PML)厚度為50個(gè)網(wǎng)格.算法采用時(shí)間上二階精度、空間上十二階精度的交錯(cuò)網(wǎng)格差分格式,空間步長為10m,時(shí)間步長為0.5ms.圖4為模擬的不同時(shí)刻的震源波場快照,震源位于模型的中間位置,時(shí)間總采樣點(diǎn)數(shù)為1000.圖4a為0.25s的正傳震源波場快照,圖4b為最大時(shí)刻0.5s的震源波場快照,圖4c為利用本文提出的改進(jìn)的PML邊界條件從最大時(shí)刻0.5s逆時(shí)回推的0.25s的波場快照.對(duì)比可知,模型內(nèi)部計(jì)算區(qū)域的波場可以通過最后時(shí)刻的波場逆時(shí)回推完整的恢復(fù)出來.2.2減少監(jiān)控存儲(chǔ)邊界波場相對(duì)于CPU,GPU在處理能力和存儲(chǔ)器帶寬上有明顯的優(yōu)勢,GPU可以通過增加并行處理單元和存儲(chǔ)器控制單元來提高處理能力和存儲(chǔ)器帶寬.GPU加速的性能受并行算法和GPU硬件構(gòu)架的影響,程序優(yōu)化程度的好壞影響計(jì)算效率.GPU通過PCI-E總線和主機(jī)相連,GPU輸入和輸出的吞吐量受到IO帶寬的限制.由于PCI-E帶寬相對(duì)較小,為提高程序的計(jì)算效率,應(yīng)該盡量減小CPU和GPU間數(shù)據(jù)傳輸量.當(dāng)需要存儲(chǔ)的邊界波場較小時(shí),可以把邊界波場全部存儲(chǔ)在顯存中,在設(shè)備端(GPU)進(jìn)行逆時(shí)回推計(jì)算時(shí),直接用存儲(chǔ)在設(shè)備端的邊界波場代替逆時(shí)回推過程中的邊界波場;當(dāng)時(shí)間離散點(diǎn)數(shù)很大時(shí),需要存儲(chǔ)的邊界波場會(huì)超出設(shè)備端(GPU)的顯存容量,需要把邊界波場存儲(chǔ)在主機(jī)端內(nèi)存中,通過PCI-E總線進(jìn)行CPU和GPU間數(shù)據(jù)傳輸.為了減小CPU和GPU間的數(shù)據(jù)傳輸量,提高計(jì)算效率,可以采取如下策略:在正演過程中設(shè)置NC個(gè)checkpoint時(shí)間點(diǎn),在每個(gè)checkpoint時(shí)間點(diǎn)處保存該時(shí)刻的全部波場,把這些checkpoint時(shí)間點(diǎn)處的波場保存到主機(jī)端內(nèi)存中,在設(shè)備端只保存相鄰兩個(gè)checkpoint時(shí)間點(diǎn)之間的邊界波場,實(shí)現(xiàn)過程如圖5.對(duì)比于把全部邊界波場保存于主機(jī)端內(nèi)存這種方法,圖5所示的這種策略可以減小CPU和GPU間的數(shù)據(jù)傳輸量,同時(shí)沒有增加計(jì)算量.針對(duì)公式(2),我們?cè)赾heckpoint離散時(shí)間點(diǎn)處,需要保存的完整波場有P、v3模型實(shí)驗(yàn)計(jì)算3.1逆時(shí)回推波場并加速正演模擬圖7a為Hess2DVTI速度模型,逆時(shí)偏移算法采用了交錯(cuò)網(wǎng)格有限差分格式,時(shí)間采用二階精度、空間采用十二階精度,網(wǎng)格間距為dx=10m、dz=10m,時(shí)間步長為0.3ms,震源位置位于地表18千米處.圖7b、圖7c為震源波場正傳1.5s時(shí)刻和采用改進(jìn)PML邊界存儲(chǔ)策略逆時(shí)回推1.5s時(shí)刻波場快照對(duì)比圖.分析圖7b和圖7c可以看出,有效PML邊界存儲(chǔ)策略可以通過最后時(shí)刻的波場逆時(shí)回推完整的恢復(fù)內(nèi)部模擬計(jì)算區(qū)域的波場.沿圖7b和7c中紅線位置,抽取快照波場的一部分對(duì)比分析,如圖7d和7e所示,改進(jìn)的PML邊界存儲(chǔ)策略,可以使模型內(nèi)部計(jì)算區(qū)域的波場完整的恢復(fù)出來,同時(shí)恢復(fù)出來的震源波場同震源正向傳播的波場振幅差可以忽略,幾乎沒有能量損失.表1為離散采樣時(shí)間點(diǎn)數(shù)nt=1000時(shí),上述Hess2DVTI模型正演模擬分別采用CPU計(jì)算和CPU/GPU并行計(jì)算的耗時(shí)對(duì)比.從表1中可以看出,CPU/GPU并行加速計(jì)算技術(shù)可以使Hess2DVTI模型的正演效率提高100倍.3.2hessdvti模型Hess2DVTI模型左部有一個(gè)高速鹽丘體,中部有有尖滅構(gòu)造,右部有較為陡峭的傾斜斷層,中部的尖滅構(gòu)造是成像難點(diǎn).模型參數(shù)包括縱波速度、各向異性參數(shù)ε、變異系數(shù)δ.逆時(shí)偏移算法采用時(shí)間二階、空間十二階的交錯(cuò)網(wǎng)格有限差分格式,網(wǎng)格間距dx和dz均為10m,時(shí)間步長為0.3ms.該模型的縱波速度如圖7a所示,圖8a為模型各向異性參數(shù)ε,圖8b為模型變異系數(shù)δ,圖8c為單炮逆時(shí)偏移成果圖、其中炮點(diǎn)位于地表18km處,圖8d為圖8c經(jīng)過拉普拉斯算子濾波去除低頻噪音的結(jié)果圖,圖8e為2DHessVTI模型360炮疊加的逆時(shí)偏移成果圖.表2為取Hess2DVTI模型的大小為3000×1500網(wǎng)格,離散時(shí)間點(diǎn)數(shù)nt=15000時(shí)幾種不同的逆時(shí)偏移存儲(chǔ)策略的計(jì)算量和存儲(chǔ)量的對(duì)比,其中NC代表checkpoint時(shí)間點(diǎn)數(shù).可以看出有效PML邊界存儲(chǔ)顯著地降低了波場的存儲(chǔ)量.4基于gpu加速技術(shù)的逆時(shí)偏移算法疊前