宋利偉,石 穎,柯 璇,王維紅,李松齡
(東北石油大學(xué),黑龍江大慶163318)
基于雙程波動(dòng)方程的逆時(shí)偏移方法,無(wú)地層傾角限制,適應(yīng)于各種波現(xiàn)象,是復(fù)雜地下構(gòu)造成像的有效技術(shù)[1-3]。由于沉積地層中砂泥巖薄互層、灰?guī)r頁(yè)巖薄互層以及大塊巖石中定向排列的垂直裂隙等常常引起地下地層各向異性[4],在各向異性介質(zhì)中地震波彈性特征隨方向發(fā)生變化,具體表現(xiàn)為地震波傳播速度隨著傳播方向的變化而變化,所以如果假設(shè)地下介質(zhì)仍為各向同性,勢(shì)必影響地震波場(chǎng)模擬精度和逆時(shí)偏移成像分辨率。雖然各向異性介質(zhì)彈性波偏移理論已經(jīng)取得了一定突破[5-7],但在應(yīng)用過(guò)程中依然會(huì)遇到障礙。其一,彈性波是由一個(gè)三分量的位移矢量表征,包括兩種模式即縱波和橫波,如利用有限差分法求解波場(chǎng),高昂的計(jì)算成本和巨大的存儲(chǔ)量令人難以接受;其二,考慮橫波速度較慢,為了滿足波場(chǎng)穩(wěn)定延拓條件,須采用較小時(shí)間步長(zhǎng)和空間網(wǎng)格間距,必然會(huì)影響計(jì)算效率;其三,縱波與橫波分離以及實(shí)際彈性參數(shù)的確定等都是其制約因素[8]。因此,對(duì)于地震波場(chǎng)在各向異性介質(zhì)傳播特征的研究聚焦于縱波[9-12]。
ALKHALIFAH[13]通過(guò)聲學(xué)近似假設(shè)提出適用于VTI介質(zhì)的四階波動(dòng)方程,該方程qP波特征與彈性波方程縱波分量特征吻合。但是,利用四階波動(dòng)方程不便于波場(chǎng)的求解,ZHOU等[14]引入一個(gè)輔助變量,將qP波方程從四階簡(jiǎn)化成兩個(gè)耦合的二階波動(dòng)方程,從而高效地對(duì)qP波進(jìn)行數(shù)值模擬?;谏鲜鲅芯?FLETCHER等[15-16]通過(guò)坐標(biāo)變換將VTI介質(zhì)qP波方程推廣為T(mén)TI介質(zhì)波動(dòng)方程。隨后,許多研究學(xué)者在ALKHALIFAH提出的四階qP波方程的基礎(chǔ)上,分別引入不同的輔助變量,得到了一系列的二階波動(dòng)方程[17-19]。李博等[20]針對(duì)TTI介質(zhì)逆時(shí)偏移算法的穩(wěn)定性進(jìn)行了詳盡的探討。
對(duì)于現(xiàn)有的二階qP波方程,盡管將方程中橫波速度設(shè)置為0,但是地震波場(chǎng)模擬仍出現(xiàn)偽橫波分量,并且,當(dāng)?shù)刭|(zhì)Thomsen參數(shù)變化劇烈時(shí),qP波方程中的偽橫波嚴(yán)重影響數(shù)值模擬的穩(wěn)定性,進(jìn)而導(dǎo)致逆時(shí)偏移方法難以得到理想的成像結(jié)果。TSVANKIN[21]提出改變震源處Thomsen參數(shù),雖然壓制了偽橫波成分,卻改變了真實(shí)地下介質(zhì)參數(shù);ZHANG[22]提出的濾波投影法壓制偽橫波取得了理想效果;CHENG等[23]從物理學(xué)的角度,利用波矢量方向與偽橫波震動(dòng)方向的偏差,從耦合的波場(chǎng)中提取出標(biāo)量的qP波;PESTANA[24]利用精確頻散關(guān)系,推導(dǎo)出波場(chǎng)的遞推關(guān)系,得到了解耦P波場(chǎng)。這些壓制偽橫波的方法提升了逆時(shí)偏移成像精度。
如果采用互相關(guān)成像條件實(shí)現(xiàn)TTI介質(zhì)逆時(shí)偏移,需要保存所有時(shí)刻的震源波場(chǎng)或檢波點(diǎn)波場(chǎng),因此,對(duì)存儲(chǔ)空間的需求較大。SYMES[25]提出checkpiont方法,僅存儲(chǔ)部分時(shí)刻的地震波場(chǎng),互相關(guān)成像時(shí)再進(jìn)行波場(chǎng)重建,但是該方法仍然需要權(quán)衡計(jì)算效率和存儲(chǔ)成本。從另一個(gè)角度講,王保利等[26]提出保存若干PML邊界層,再進(jìn)行波場(chǎng)重建的存儲(chǔ)策略,并將該方法應(yīng)用到了聲波逆時(shí)偏移中,但是該方法在重建時(shí)利用了PML邊界層中經(jīng)過(guò)衰減的波場(chǎng)信息,理論上存在重建誤差。
在計(jì)算機(jī)技術(shù)取得長(zhǎng)足進(jìn)步的背景下,GPU技術(shù)為高效的科學(xué)計(jì)算提供了強(qiáng)有力的工具。為了提高逆時(shí)偏移的計(jì)算效率,逆時(shí)偏移算法與GPU硬件相結(jié)合成為了歷史的必然,形成的GPU/CPU并行加速計(jì)算方案滿足了工業(yè)化生產(chǎn)的需求[27-29]。
本文在前人研究的基礎(chǔ)上,考慮地下介質(zhì)的各向異性,研究并實(shí)現(xiàn)了TTI介質(zhì)波場(chǎng)傳播特征的精確描述及逆時(shí)偏移成像方法,并用BP2007模型數(shù)據(jù)進(jìn)行了測(cè)試。
刻畫(huà)地震波在各向異性介質(zhì)中傳播的3個(gè)要素是波的震動(dòng)方向、群速度和相速度。對(duì)于二維VTI介質(zhì),TSVANKIN[21]通過(guò)求解Christoffel方程得到qP-qSV波的精準(zhǔn)相速度公式,即:
(1a)
(1b)
式中:vP(φ)是P波相速度;φ是相角;vPz和vSz分別為縱波和橫波沿著對(duì)稱軸方向的傳播速度;ε和δ為T(mén)homsen參數(shù);“±”中正號(hào)對(duì)應(yīng)qP波,負(fù)號(hào)對(duì)應(yīng)qSV波。利用相速度、頻率和波數(shù)之間的關(guān)系,可以得到VTI介質(zhì)色散方程,形式如下:
式中:ω為角頻率;kx和kz為波矢量;vPx為縱波沿垂直對(duì)稱軸方向的速度;vPn為動(dòng)校正速度。假設(shè)TTI介質(zhì)的對(duì)稱軸與垂直方向夾角為θ,(2)式中波矢量通過(guò)坐標(biāo)變換為:
(3)
將坐標(biāo)變換后的波矢量代入(2)式有:
(5)
(6)
對(duì)方程(5)和方程(6)進(jìn)行反傅里葉變換就將四階微分方程轉(zhuǎn)化為適用于TTI介質(zhì)的兩個(gè)二階耦合微分方程:
對(duì)于各向同性介質(zhì),方程(7)中的兩式等價(jià),其二階微分算子形式如下:
(8a)
(8b)
如果將(8)式中θ設(shè)置為0,那么TTI介質(zhì)方程就退化為VTI介質(zhì)方程。
如果利用互相關(guān)成像條件,TTI介質(zhì)逆時(shí)偏移成像過(guò)程需要對(duì)每一時(shí)刻的震源波場(chǎng)和檢波點(diǎn)波場(chǎng)做互相關(guān)。因此,地震波場(chǎng)的模擬精度直接影響最后成像結(jié)果的分辨率,由于數(shù)值模擬只能在有限區(qū)域內(nèi)計(jì)算,存在截?cái)噙吔?而邊界的內(nèi)外兩側(cè)具有不同的波阻抗,因而相當(dāng)于人為增加了一個(gè)反射界面,從而在界面產(chǎn)生虛假反射,導(dǎo)致結(jié)果錯(cuò)誤[30-31]。因此,需要為偏移區(qū)域添加吸收層以避免能量反射。本文采用經(jīng)典PML邊界[32]吸收到達(dá)邊界的波場(chǎng)能量。根據(jù)聲波PML吸收邊界條件的架構(gòu)[33-35],將(7)式變換為下列形式:
(9a)
(9b)
式中:d(x,z)為衰減因子。為了能夠?qū)崿F(xiàn)數(shù)值模擬,時(shí)間方向上的二階和一階導(dǎo)數(shù)項(xiàng)采用公式(10)和公式(11)進(jìn)行計(jì)算。
將(10)式和(11)式代入(9)式,得到TTI介質(zhì)qP波方程有限差分格式如下:
由公式(12)可以得到波場(chǎng)重建的逆向延拓表達(dá)式,如公式(13)所示。
采用高階中心有限差分方法求解(8)式中的空間導(dǎo)數(shù)項(xiàng),水平方向2N階差分,垂直方向2M階差分,形式如下:
(14d)
式中:an,bm,cn,dm為差分系數(shù)。
如果采用聲學(xué)近似的方法,將方程(12)中橫波沿對(duì)稱軸的速度設(shè)置為0,雖然能夠提高計(jì)算效率,但因橫波沿著其它方向的速度并不為0,仍有低振幅低波速的qSV波成分。當(dāng)TTI介質(zhì)的對(duì)稱軸傾角較大或地下介質(zhì)的對(duì)稱軸傾角變化較大時(shí),qSV波的存在將導(dǎo)致波場(chǎng)傳播不穩(wěn)定。雖然通過(guò)修改各向異性參數(shù)能夠消除頻散的影響,但也會(huì)改變波場(chǎng)傳播的動(dòng)力學(xué)特征。為了避免地震波場(chǎng)異常傳播,并且能夠準(zhǔn)確地描述波場(chǎng)運(yùn)動(dòng)學(xué)傳播特征,在TTI介質(zhì)波動(dòng)方程中引入σ參數(shù),其數(shù)學(xué)描述如下:
(15)
σ參數(shù)在很大程度上決定了qSV波的運(yùn)動(dòng)學(xué)特征。通過(guò)數(shù)值模擬,當(dāng)設(shè)置σ=0.75時(shí),頻散得以壓制,確保波場(chǎng)的穩(wěn)定傳播。
基于互相關(guān)成像條件的逆時(shí)偏移方法需要存儲(chǔ)全部時(shí)刻震源波場(chǎng),或者存儲(chǔ)全部時(shí)刻的檢波點(diǎn)波場(chǎng),因而存儲(chǔ)成本巨大。以BP2007模型為例,重建試算截取模型水平方向2000個(gè)網(wǎng)格點(diǎn),垂直方向1801個(gè)網(wǎng)格點(diǎn),地震記錄總時(shí)長(zhǎng)9.2s,為了滿足波場(chǎng)延拓穩(wěn)定性條件,選取時(shí)間步長(zhǎng)0.5ms,需要延拓18401次,保存所有震源波場(chǎng)需要內(nèi)存約246.9GB。因此,SHI等[36]考慮存儲(chǔ)較少的波場(chǎng)信息,進(jìn)而重建震源波場(chǎng);王保利等[26]利用存儲(chǔ)若干PML吸收層信息實(shí)現(xiàn)了波場(chǎng)重建并將該方法應(yīng)用到聲波逆時(shí)偏移成像中。對(duì)于TTI介質(zhì),逆時(shí)偏移成
像對(duì)數(shù)據(jù)存儲(chǔ)需求更高,為了增強(qiáng)方法的實(shí)用性,本文中提出了實(shí)用的TTI介質(zhì)的checkpoint波場(chǎng)重建方法。
假定地震波場(chǎng)水平方向?yàn)镹x個(gè)網(wǎng)格,垂直方向?yàn)镹z個(gè)網(wǎng)格,為了避免邊界反射噪聲,在有效波場(chǎng)外側(cè)添加PML吸收衰減層(圖1中綠色區(qū)域),在有效波場(chǎng)內(nèi)部設(shè)置內(nèi)邊界層(圖1中藍(lán)色區(qū)域)。內(nèi)邊界層數(shù)設(shè)定為有限差分階數(shù)的一半,如果采用空間12階差分,則內(nèi)邊界為6層。各向異性介質(zhì)checkpoint波場(chǎng)重建由兩次震源波場(chǎng)正傳和一次波場(chǎng)反傳組成。其中一次震源波場(chǎng)正傳在最小時(shí)刻和最大時(shí)刻之間進(jìn)行,而另一次震源波場(chǎng)正傳和一次震源波場(chǎng)反傳逐個(gè)在兩個(gè)checkpoint之間進(jìn)行。具體如下。
圖1 波場(chǎng)存儲(chǔ)示意
假設(shè)震源波場(chǎng)正向傳播的時(shí)間點(diǎn)為N×m,N為在此傳播過(guò)程中確定的checkpoint個(gè)數(shù),m為每?jī)蓚€(gè)checkpoint間的時(shí)間點(diǎn)數(shù)。在波場(chǎng)正向延拓過(guò)程中,保存每個(gè)checkpoint及其前一時(shí)間點(diǎn)處的p和q全部波場(chǎng)(圖1中PML吸收層、內(nèi)邊界層和內(nèi)部區(qū))至計(jì)算機(jī)內(nèi)存,圖2中黑點(diǎn)表示存儲(chǔ)p和q全部波場(chǎng)的時(shí)間位置。例如,當(dāng)波場(chǎng)延拓至checkpointN時(shí),記錄(N×m-1)和N×m時(shí)間點(diǎn)對(duì)應(yīng)的p和q全部波場(chǎng)。
以checkpoint(N-1)至checkpointN期間的波場(chǎng)重建為例,需從計(jì)算機(jī)內(nèi)存中讀取(N-1)×m-1及(N-1)×m時(shí)間點(diǎn)對(duì)應(yīng)的p和q全部波場(chǎng)信息作為初值進(jìn)行m次正向延拓,并將每次延拓波場(chǎng)的內(nèi)邊界層(圖1中藍(lán)色區(qū)域)信息保存至GPU的顯存。利用公式(13)將N×m-1和N×m節(jié)點(diǎn)對(duì)應(yīng)的p和q全部波場(chǎng)作為反向延拓的初值,每延拓一次將上一步驟中保存至顯存上p和q的內(nèi)邊界層(圖1中藍(lán)色區(qū)域)信息載入到重建波場(chǎng)的對(duì)應(yīng)位置,延拓m次便可重建出checkpoint(N-1)至checkpointN之間對(duì)應(yīng)的震源波場(chǎng)。其它任兩個(gè)checkpoint之間波場(chǎng)重建及內(nèi)邊界層的保存辦法同上,進(jìn)而重建出N×m至最小時(shí)刻的震源波場(chǎng)。本文提出的TTI介質(zhì)逆時(shí)偏移成像方法是在checkpoint框架下,通過(guò)保存內(nèi)邊界層信息實(shí)現(xiàn)震源波場(chǎng)重建,而非保存經(jīng)過(guò)吸收衰減后的PML吸收層波場(chǎng)信息,因此確保了波場(chǎng)重建精度。另外,在波場(chǎng)重建過(guò)程中,與保存所有波場(chǎng)相比,僅保存兩個(gè)checkpoint之間波場(chǎng)的內(nèi)邊界層信息,雖然多了一次波場(chǎng)正向延拓的計(jì)算量,但是GPU/CPU協(xié)同并行計(jì)算技術(shù)在本文方法中的應(yīng)用,較大地提高了計(jì)算效率,充分補(bǔ)償了由此帶來(lái)的計(jì)算負(fù)擔(dān),達(dá)到了以計(jì)算換存儲(chǔ)的目的。需要指出的是在延拓時(shí)間點(diǎn)一定的情況下,checkpoint越多,GPU顯存占用就越少,雖然降低了對(duì)GPU硬件的要求,但是增加了計(jì)算機(jī)內(nèi)存與顯存之間數(shù)據(jù)交換的次數(shù),這將導(dǎo)致計(jì)算耗時(shí)明顯增長(zhǎng),所以在確保兩checkpoint間的波場(chǎng)內(nèi)邊界占用空間不超過(guò)GPU顯存的情況下,可適當(dāng)減少checkpoint數(shù),以降低數(shù)據(jù)交換頻率,從而達(dá)到計(jì)算效率的優(yōu)化。
圖2 checkpoint存儲(chǔ)示意
為了驗(yàn)證上述推導(dǎo)帶有PML吸收邊界qP波方程有限差分格式的正確有效性,測(cè)試了均勻TTI介質(zhì)的地震波場(chǎng)傳播。假定模型在水平和垂直方向的網(wǎng)格數(shù)均為320,水平方向和垂直方向的網(wǎng)格間距均為6.25m,沿對(duì)稱軸方向的縱波速度vPz=3000m/s,Thomsen各向異性參數(shù)δ=0.1,ε=0.3,PML吸收邊界層數(shù)為50,震源位于模型的中央,主頻為30Hz的雷克子波作為激發(fā)震源,時(shí)間采樣間隔為0.5ms。
利用本文推導(dǎo)的公式(12)進(jìn)行數(shù)值模擬得到的波場(chǎng)快照如圖3所示。圖3a為對(duì)稱軸傾角θ=0,沿著對(duì)稱軸方向的橫波速度為0,波傳播時(shí)間為200ms時(shí)p的波場(chǎng)快照。圖3b為對(duì)稱軸傾角θ=0,沿著對(duì)稱軸方向的橫波速度為0,波傳播時(shí)間為350ms時(shí)p的波場(chǎng)快照,觀察圖3b可知,PML邊界對(duì)地震邊界波場(chǎng)的吸收效果較為理想。圖3c為對(duì)稱軸傾角θ=30°,沿著對(duì)稱軸方向的橫波速度為0,波傳播時(shí)間為200ms時(shí)p的波場(chǎng)快照。從圖3c可以看出,qSV波附近位置出現(xiàn)嚴(yán)重頻散現(xiàn)象。圖3d為對(duì)稱軸傾角θ=30°,沿著對(duì)稱軸方向的橫波速度為1500m/s,波傳播時(shí)間為200ms時(shí)的p的波場(chǎng)快照,可見(jiàn),加入適當(dāng)?shù)臋M波速度或者引入σ參數(shù)后,頻散得到了有效抑制,TTI介質(zhì)數(shù)值模擬精度較高。圖3a,圖3b和圖3c 中雖然橫波沿對(duì)稱軸方向的速度為0,但是沿著其它方向的橫波速度并不為0,因此出現(xiàn)了qSV波,如圖中菱形狀波場(chǎng)所示。
圖3 p的波場(chǎng)快照a vSz=0,θ=0,t=200ms; b vSz=0,θ=0,t=350ms; c vSz=0,θ=30°,t=200ms; d vSz=1500m/s,θ=30°,t=200ms
測(cè)試結(jié)果表明,對(duì)于圖3d所示的數(shù)值模擬,采用CPU算法共耗時(shí)7min,采用CPU/GPU協(xié)同并行加速算法耗時(shí)8.5s,計(jì)算效率提高約50倍。如果將模型的橫、縱向網(wǎng)格點(diǎn)數(shù)由320增加至3200,數(shù)據(jù)量增加100倍,采用CPU/GPU協(xié)同并行加速算法耗時(shí)40.9s,計(jì)算耗時(shí)是原數(shù)據(jù)的5倍。由此可見(jiàn),在波場(chǎng)正演模擬中引入CPU/GPU協(xié)同并行加速技術(shù),可以大幅提升計(jì)算效率,隨著模型網(wǎng)格點(diǎn)數(shù)的增大,并行計(jì)算耗時(shí)也相應(yīng)地增加,但并不是線性增加。
利用BP2007復(fù)雜構(gòu)造模型檢驗(yàn)本文介紹的TTI介質(zhì)逆時(shí)偏移成像方法。模型水平方向12596個(gè)網(wǎng)格點(diǎn),垂直方向1801網(wǎng)格點(diǎn),網(wǎng)格間距6.25m,水平長(zhǎng)度78.70km,垂直深度11.25km。該模型數(shù)據(jù)共1641炮,每炮800道,最小偏移距37.5m,道間距12.5m,時(shí)間方向1151個(gè)采樣點(diǎn),采樣間隔為8ms。本文截取其水平方向43.70km,共626炮數(shù)據(jù)進(jìn)行TTI介質(zhì)逆時(shí)偏移成像測(cè)試,截取部分模型參數(shù)如圖4所示。
首先利用截取部分模型來(lái)測(cè)試基于checkpoint方法重建qP波的可行性和有效性。截取模型的水平方向2000個(gè)網(wǎng)格點(diǎn),垂直方向1801網(wǎng)格點(diǎn),時(shí)間采樣間隔為0.5ms,利用主頻為30Hz的雷克子波作為激發(fā)震源,激發(fā)點(diǎn)位于水平方向1000網(wǎng)格點(diǎn)、垂直方向900網(wǎng)格點(diǎn)處。震源波場(chǎng)傳播至1.5s時(shí)的波場(chǎng)快照如圖5a所示,震源波場(chǎng)傳播至最大時(shí)刻后,再根據(jù)本文提出的checkpoint方法和內(nèi)邊界保存技術(shù)對(duì)震源波場(chǎng)進(jìn)行分時(shí)段的正向傳播和逆向傳播,逆向傳播至1.5s時(shí)的震源波場(chǎng)如圖5b所示。正向傳播的震源波場(chǎng)和逆向傳播的震源重建波場(chǎng)在1.5s時(shí)刻的波場(chǎng)誤差如圖5c所示,誤差波場(chǎng)振幅與有效波場(chǎng)振幅相差5個(gè)數(shù)量級(jí),可忽略不計(jì)。由此可知,在checkpoint方法框架下保存波場(chǎng)內(nèi)邊界的方法,可以有效地進(jìn)行波場(chǎng)重建,豐富了逆時(shí)偏移算法的研究,為T(mén)TI介質(zhì)逆時(shí)偏移方法的工業(yè)化應(yīng)用帶來(lái)了更廣闊的發(fā)展空間。
利用截取的BP2007模型進(jìn)行偏移試算。單炮偏移區(qū)域水平方向2000個(gè)網(wǎng)格點(diǎn),垂直方向1801個(gè)網(wǎng)格點(diǎn),空間網(wǎng)格間距為6.25m,時(shí)間采樣間隔為0.5ms。先將震源子波和地層模型對(duì)震源波場(chǎng)正向延拓,保存震源波場(chǎng)部分信息,然后根據(jù)保存的波場(chǎng)信息,基于checkpoint方法對(duì)震源波場(chǎng)重建,與此同時(shí)將炮集數(shù)據(jù)加載到檢波點(diǎn)得到沿時(shí)間逆向的檢波點(diǎn)波場(chǎng),最后經(jīng)震源歸一化成像條件得到單炮偏移剖面。經(jīng)單炮逆時(shí)偏移算法測(cè)試,采用CPU/GPU協(xié)同并行加速算法耗時(shí)39min,從上述波場(chǎng)重建方法可知該時(shí)長(zhǎng)受checkpoint數(shù)的影響,而采用CPU算法耗時(shí)1755min,GPU硬件的應(yīng)用可大幅提升逆時(shí)偏移的計(jì)算效率。對(duì)單炮偏移剖面利用拉普拉斯方法去除低頻噪聲后,將626炮偏移結(jié)果進(jìn)行多炮疊加得到最后的偏移剖面,如圖6所示。從圖6可知,基于qP波方程得到的剖面同相軸連續(xù)性好,兩個(gè)陡傾角鹽丘與模型位置匹配。
圖4 TTI介質(zhì)模型參數(shù)a 縱波沿著對(duì)稱軸的速度; b Thomsen參數(shù)θ; c Thomsen參數(shù)δ; d Thomsen參數(shù)ε
為了便于分析本文逆時(shí)偏移方法的精確性,根據(jù)偏移速度場(chǎng),結(jié)合常密度參數(shù),計(jì)算得到反射系數(shù)模型,如圖7所示。分別在圖6和圖7的水平方向11.25km位置抽取單道數(shù)據(jù)如圖8所示(圖中橫軸為剖面的垂直深度,縱軸是經(jīng)歸一化的單道振幅信息)。由圖8可知,逆時(shí)偏移剖面包含了強(qiáng)反射系數(shù)信息,同時(shí)弱反射系數(shù)也得以體現(xiàn)。因此,本文的逆時(shí)偏移方法可以實(shí)現(xiàn)對(duì)地下各向異性構(gòu)造體的精細(xì)刻畫(huà)。
圖5 震源波場(chǎng)重建a 正傳震源波場(chǎng); b 重建震源波場(chǎng); c 重建震源波場(chǎng)誤差
圖6 逆時(shí)偏移剖面
圖7 反射系數(shù)模型
圖8 水平方向11.25km處單道數(shù)據(jù)對(duì)比
從精準(zhǔn)相速度公式出發(fā),推導(dǎo)了基于有限差分?jǐn)?shù)值解法的TTI介質(zhì)qP波方程的PML吸收邊界格式。由波場(chǎng)快照可知,當(dāng)波場(chǎng)能量傳播至邊界層后吸收效果良好。在波動(dòng)方程中引入橫波速度,數(shù)值試驗(yàn)結(jié)果表明,該方法有效地壓制了qSV波頻散,能夠解決波場(chǎng)傳播不穩(wěn)定的問(wèn)題。為了解決算法無(wú)法在有限顯存GPU端運(yùn)行的問(wèn)題,本文在各向異性介質(zhì)偏移成像過(guò)程中應(yīng)用了基于checkpoint方法和保存內(nèi)邊界波場(chǎng)的存儲(chǔ)策略,極大地降低了存儲(chǔ),雖然該策略增加了一次波場(chǎng)正向延拓計(jì)算量以及硬件之間數(shù)據(jù)的傳輸時(shí)間,但是偏移成像過(guò)程中應(yīng)用了GPU/CPU并行加速技術(shù),總體上大幅度提升了計(jì)算效率,實(shí)現(xiàn)以計(jì)算換存儲(chǔ)的目的,從而降低對(duì)計(jì)算硬件的依賴,為T(mén)TI介質(zhì)逆時(shí)偏移成像的工業(yè)化應(yīng)用開(kāi)辟了有效的思路。
[1]BAYSAL E,KOSLOFF D,SHERWOOD J.Reverse time migration[J].Geophysics,1983,48(11):1514-1524
[2]MCMECHAN G A.Migration by extrapolation of time-dependent boundary values[J].Geophysical Prospecting,1983,31(3):413-420
[3]WHITMORE N D.Iterative depth migration by backward time propagation[J].Expanded Abstracts of 53rdAnnual Internat SEG Mtg,1983:382-385
[4]THOMSEN L.Weak elastic anisotropy[J].Geophysics,1986,51(10):1954-1966
[5]YAN J,SAVA P.3D elastic wave mode separation for TTI media[J].Expanded Abstracts of 79thAnnual Internat SEG Mtg,2009:1197-1201
[6]YAN J,SAVA P.Elastic wave-mode separation for VTI media[J].Geophysics,2009,74(5):1954-1966
[7]HU X L,JIA X F,ZHANG W.A dynamic lattice method for elastic wave modeling and migration in TTI media[J].Expanded Abstracts of 86thAnnual Internat SEG Mtg,2016:447-451
[8]秦海旭,吳國(guó)忱.TTI介質(zhì)彈性波隨機(jī)邊界逆時(shí)偏移的實(shí)現(xiàn)[J].石油物探,2014,53(5):570-578
QIN H X,WU G C.The implementation of elastic reverse time migration in TTI media based on random boundary[J].Geophysical Prospecting for Petroleum,2014,53(5):570-578
[9]張千祥,王德利,周進(jìn)舉.二維TTI介質(zhì)的純P波波動(dòng)方程數(shù)值模擬[J].石油物探,2015,54(5):485-492
ZHANG Q X,WANG D L,ZHOU J J.Acoustic wave equation numerical simulation for pure P-wave in 2D TTI medium[J].Geophysical Prospecting for Petroleum,2015,54(5):485-492
[10]張衡,劉洪,李博,等.TTI介質(zhì)聲波方程分裂式PML吸收邊界條件研究[J].石油物探,2017,56(3):349-361
ZHANG H,LIU H,LI B,et al.The research on split PML absorbing boundary conditions of acoustic equation for TTI media[J].Geophysical Prospecting for Petroleum,2017,56(3):349-361
[11]ALKHALIFAH T.Acoustic approximations media for processing in transversely isotropic[J].Geophysics,1998,63(2):623-631
[12]HAN Q,WU R S.One-way dual-domain propagator for scalar qP-waves in VTI media[J].Geophysics,2005,70(2):D9-D17
[13]ALKHALIFAH T.An acoustic wave equation for anisotropic media[J].Geophysics,2000,65(4):1239-1250
[14]ZHOU H B,ZHANG G Q,BLOOR R.An anisotropic acoustic wave equation for modeling and migration in 2D TTI media[J].Expanded Abstracts of 76thAnnual Internat SEG Mtg,2006:194-198
[15]FLETCHER R P,DU X,FOWLER P J.Reverse time migration in tilted transversely isotropic(TTI) media[J].Geophysics,2009,74(6):WCA179-WCA187
[16]FLETCHER R P,DU X,FOWLER P J.Stabilizing acoustic reverse-time migration in TTI media[J].Expanded Abstracts of 79thAnnual Internat SEG Mtg,2009:2985-2989
[17]DUVENECK E,MILCIK P,BAKKER P M.Acoustic VTI wave equations and their application for anisotropic reverse-time migration[J].Expanded Abstracts of 78thAnnual Internat SEG Mtg,2008:2186-2190
[18]FOWLER P J,DU X,FLETCHER R P.Coupled equations for reverse-time migration in transversely isotropic media[J].Geophysics,2010,75(1):S11-S22
[19]ZHANG L B,RECTOR J W,HOVERSTEN G M.An acoustic wave equation for modeling in tilted TI media[J].Expanded Abstracts of 73rdAnnual Internat SEG Mtg,2003:153-156
[20]李博,李敏,劉紅偉,等.TTI介質(zhì)有限差分逆時(shí)偏移的穩(wěn)定性探討[J].地球物理學(xué)報(bào),2014,57(4):1366-1375
LI B,LI M,LIU H W,et al.Stability of reverse time migration in TTI media[J].Chinese Journal of Geophysics,2014,57(4):1366-1375
[21]TSVANKIN I.Seismic signatures and analysis of reflection data in anisotropic media[M].New York:Elsevier,2001:1-454
[22]ZHANG H Z.Removing S-wave noise in TTI reverse time migration[J].Expanded Abstracts of 79thAnnual Internat SEG Mtg,2009:2849-2853
[23]CHENG J B,KANG W.Simulating propagation of separated wave modes in general anisotropic media,part Ⅱ:qS-wave propagators[J].Geophysics,2016,81(2):C39-C52
[24]PESTANA R C.Separate P-and SV-wave equations for VTI media[J].12thInternational Congress of the Brazilian Geophysical Society & EXPOGEF SEG,2011:1227-1232
[25]SYMES W W.Reverse time migration with optimal checkpointing[J].Geophysics,2007,72(5):SM213-SM221
[26]王保利,高靜懷,陳文超,等.地震疊前逆時(shí)偏移的有效邊界存儲(chǔ)策略[J].地球物理學(xué)報(bào),2012,55(7):2412-2421
WANG B L,GAO J H,CHEN W C,et al.Efficient boundary storage strategies for seismic reverse migration[J].Chinese Journal of Geophysics,2012,55(7):2412-2421
[27]劉國(guó)峰,劉洪,李博,等.山地地震資料疊前時(shí)間偏移方法及其GPU實(shí)現(xiàn)[J].地球物理學(xué)報(bào),2009,52(12):3101-3108
LIU G F,LIU H,LI B,et al.Method of prestack time migration of seismic data of mountainous regions and its GPU implementation[J].Chinese Journal of Geophysics,2009,52(12):3101-3108
[28]李博,劉紅偉,劉國(guó)峰,等.地震疊前逆時(shí)偏移算法的CPU/GPU實(shí)施對(duì)策[J].地球物理學(xué)報(bào),2010,53(12):2938-2943
LI B,LIU H W,LIU G F,et al.Computational strategy of seismic pre-stack reverse time migration on CPU/GPU[J].Chinese Journal of Geophysics,2010,53(12):2938-2943
[29]劉紅偉,李博,劉洪,等.地震疊前逆時(shí)偏移高階有限差分算法及GPU實(shí)現(xiàn)[J].地球物理學(xué)報(bào),2010,53(7):1725-1733
LIU H W,LI B,LIU H,et al.The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation[J].Chinese Journal of Geophysics,2010,53(7):1725-1733
[30]柯璇,石穎,宋利偉,等.基于褶積完全匹配吸收邊界的聲波方程數(shù)值模擬[J].石油物探,2017,56(5):637-643
KE X,SHI Y,SONG L W,et al.Numerical modeling of acoustic wave equations based on convolutional perfectly matched layer absorbing boundary condition[J].Geophysical Prospecting for Petroleum,2017,56(5):637-643
[31]柯璇,石穎,張瑩瑩,等.地震疊前逆時(shí)偏移衰減隨機(jī)邊界條件研究[J].石油物探,2017,56(4):523-533
KE X,SHI Y,ZHANG Y Y,et al.A damped random boundary condition for prestack reverse time migration[J].Geophysical Prospecting for Petroleum,2017,56(4):523-533
[32]BERENGER J P.A Perfectly matched layer for the absorption of electromagnetic waves[J].Journal of Computational Physics,1994,114(2):185-200
[33]孫林潔,印興耀.基于PML邊界條件的高倍可變網(wǎng)格有限差分?jǐn)?shù)值模擬方法[J].地球物理學(xué)報(bào),2011,54(6):1614-1623
SUN L J,YIN X Y.A finite-difference scheme based on PML boundary condition with high power grid step variation[J].Chinese Journal of Geophysics,2011,54(6):1614-1623
[34]LIU Q H,TAO J P.The perfectly matched layer for acoustic waves in absorptive media[J].Journal of the Acoustical Society of America,1997,102(4):2072-2082
[35]CHEW W C,LIU Q H.Perfectly matched layers for elastodynamics:a new absorbing boundary condition[J].Journal of Computational Acoustics,1996,4(4):341-359
[36]SHI Y,WANG Y H.Reverse time migration of 3D vertical seismic profile data[J].Geophysics,2016,81(1):S31-S38