宋利偉,石　穎,柯　璇,王維紅,李松齡

(東北石油大學,黑龍江大慶163318)

基于雙程波動方程的逆時偏移方法,無地層傾角限制,適應于各種波現(xiàn)象,是復雜地下構(gòu)造成像的有效技術(shù)[1-3]。由于沉積地層中砂泥巖薄互層、灰?guī)r頁巖薄互層以及大塊巖石中定向排列的垂直裂隙等常常引起地下地層各向異性[4],在各向異性介質(zhì)中地震波彈性特征隨方向發(fā)生變化,具體表現(xiàn)為地震波傳播速度隨著傳播方向的變化而變化,所以如果假設地下介質(zhì)仍為各向同性,勢必影響地震波場模擬精度和逆時偏移成像分辨率。雖然各向異性介質(zhì)彈性波偏移理論已經(jīng)取得了一定突破[5-7],但在應用過程中依然會遇到障礙。其一,彈性波是由一個三分量的位移矢量表征,包括兩種模式即縱波和橫波,如利用有限差分法求解波場,高昂的計算成本和巨大的存儲量令人難以接受;其二,考慮橫波速度較慢,為了滿足波場穩(wěn)定延拓條件,須采用較小時間步長和空間網(wǎng)格間距,必然會影響計算效率;其三,縱波與橫波分離以及實際彈性參數(shù)的確定等都是其制約因素[8]。因此,對于地震波場在各向異性介質(zhì)傳播特征的研究聚焦于縱波[9-12]。

ALKHALIFAH[13]通過聲學近似假設提出適用于VTI介質(zhì)的四階波動方程,該方程qP波特征與彈性波方程縱波分量特征吻合。但是,利用四階波動方程不便于波場的求解,ZHOU等[14]引入一個輔助變量,將qP波方程從四階簡化成兩個耦合的二階波動方程,從而高效地對qP波進行數(shù)值模擬?；谏鲜鲅芯?FLETCHER等[15-16]通過坐標變換將VTI介質(zhì)qP波方程推廣為TTI介質(zhì)波動方程。隨后,許多研究學者在ALKHALIFAH提出的四階qP波方程的基礎上,分別引入不同的輔助變量,得到了一系列的二階波動方程[17-19]。李博等[20]針對TTI介質(zhì)逆時偏移算法的穩(wěn)定性進行了詳盡的探討。

對于現(xiàn)有的二階qP波方程,盡管將方程中橫波速度設置為0,但是地震波場模擬仍出現(xiàn)偽橫波分量,并且,當?shù)刭|(zhì)Thomsen參數(shù)變化劇烈時,qP波方程中的偽橫波嚴重影響數(shù)值模擬的穩(wěn)定性,進而導致逆時偏移方法難以得到理想的成像結(jié)果。TSVANKIN[21]提出改變震源處Thomsen參數(shù),雖然壓制了偽橫波成分,卻改變了真實地下介質(zhì)參數(shù);ZHANG[22]提出的濾波投影法壓制偽橫波取得了理想效果;CHENG等[23]從物理學的角度,利用波矢量方向與偽橫波震動方向的偏差,從耦合的波場中提取出標量的qP波;PESTANA[24]利用精確頻散關(guān)系,推導出波場的遞推關(guān)系,得到了解耦P波場。這些壓制偽橫波的方法提升了逆時偏移成像精度。

如果采用互相關(guān)成像條件實現(xiàn)TTI介質(zhì)逆時偏移,需要保存所有時刻的震源波場或檢波點波場,因此,對存儲空間的需求較大。SYMES[25]提出checkpiont方法,僅存儲部分時刻的地震波場,互相關(guān)成像時再進行波場重建,但是該方法仍然需要權(quán)衡計算效率和存儲成本。從另一個角度講,王保利等[26]提出保存若干PML邊界層,再進行波場重建的存儲策略,并將該方法應用到了聲波逆時偏移中,但是該方法在重建時利用了PML邊界層中經(jīng)過衰減的波場信息,理論上存在重建誤差。

在計算機技術(shù)取得長足進步的背景下,GPU技術(shù)為高效的科學計算提供了強有力的工具。為了提高逆時偏移的計算效率,逆時偏移算法與GPU硬件相結(jié)合成為了歷史的必然,形成的GPU/CPU并行加速計算方案滿足了工業(yè)化生產(chǎn)的需求[27-29]。

本文在前人研究的基礎上,考慮地下介質(zhì)的各向異性,研究并實現(xiàn)了TTI介質(zhì)波場傳播特征的精確描述及逆時偏移成像方法,并用BP2007模型數(shù)據(jù)進行了測試。

1　TTI介質(zhì)qP波方程完全匹配吸收邊界有限差分格式的推導

刻畫地震波在各向異性介質(zhì)中傳播的3個要素是波的震動方向、群速度和相速度。對于二維VTI介質(zhì),TSVANKIN[21]通過求解Christoffel方程得到qP-qSV波的精準相速度公式,即:

(1a)

(1b)

式中:vP(φ)是P波相速度;φ是相角;vPz和vSz分別為縱波和橫波沿著對稱軸方向的傳播速度;ε和δ為Thomsen參數(shù);“±”中正號對應qP波,負號對應qSV波。利用相速度、頻率和波數(shù)之間的關(guān)系,可以得到VTI介質(zhì)色散方程,形式如下:

式中:ω為角頻率;kx和kz為波矢量;vPx為縱波沿垂直對稱軸方向的速度;vPn為動校正速度。假設TTI介質(zhì)的對稱軸與垂直方向夾角為θ,(2)式中波矢量通過坐標變換為:

(3)

將坐標變換后的波矢量代入(2)式有:

(5)

(6)

對方程(5)和方程(6)進行反傅里葉變換就將四階微分方程轉(zhuǎn)化為適用于TTI介質(zhì)的兩個二階耦合微分方程:

對于各向同性介質(zhì),方程(7)中的兩式等價,其二階微分算子形式如下:

(8a)

(8b)

如果將(8)式中θ設置為0,那么TTI介質(zhì)方程就退化為VTI介質(zhì)方程。

如果利用互相關(guān)成像條件,TTI介質(zhì)逆時偏移成像過程需要對每一時刻的震源波場和檢波點波場做互相關(guān)。因此,地震波場的模擬精度直接影響最后成像結(jié)果的分辨率,由于數(shù)值模擬只能在有限區(qū)域內(nèi)計算,存在截斷邊界,而邊界的內(nèi)外兩側(cè)具有不同的波阻抗,因而相當于人為增加了一個反射界面,從而在界面產(chǎn)生虛假反射,導致結(jié)果錯誤[30-31]。因此,需要為偏移區(qū)域添加吸收層以避免能量反射。本文采用經(jīng)典PML邊界[32]吸收到達邊界的波場能量。根據(jù)聲波PML吸收邊界條件的架構(gòu)[33-35],將(7)式變換為下列形式:

(9a)

(9b)

式中:d(x,z)為衰減因子。為了能夠?qū)崿F(xiàn)數(shù)值模擬,時間方向上的二階和一階導數(shù)項采用公式(10)和公式(11)進行計算。

將(10)式和(11)式代入(9)式，得到TTI介質(zhì)qP波方程有限差分格式如下:

由公式(12)可以得到波場重建的逆向延拓表達式，如公式(13)所示。

采用高階中心有限差分方法求解(8)式中的空間導數(shù)項,水平方向2N階差分,垂直方向2M階差分,形式如下:

(14d)

式中:an,bm,cn,dm為差分系數(shù)。

如果采用聲學近似的方法,將方程(12)中橫波沿對稱軸的速度設置為0,雖然能夠提高計算效率,但因橫波沿著其它方向的速度并不為0,仍有低振幅低波速的qSV波成分。當TTI介質(zhì)的對稱軸傾角較大或地下介質(zhì)的對稱軸傾角變化較大時,qSV波的存在將導致波場傳播不穩(wěn)定。雖然通過修改各向異性參數(shù)能夠消除頻散的影響,但也會改變波場傳播的動力學特征。為了避免地震波場異常傳播,并且能夠準確地描述波場運動學傳播特征,在TTI介質(zhì)波動方程中引入σ參數(shù),其數(shù)學描述如下:

(15)

σ參數(shù)在很大程度上決定了qSV波的運動學特征。通過數(shù)值模擬,當設置σ=0.75時,頻散得以壓制,確保波場的穩(wěn)定傳播。

2　基于checkpoint方法的TTI介質(zhì)逆時偏移波場重建策略

基于互相關(guān)成像條件的逆時偏移方法需要存儲全部時刻震源波場,或者存儲全部時刻的檢波點波場,因而存儲成本巨大。以BP2007模型為例,重建試算截取模型水平方向2000個網(wǎng)格點,垂直方向1801個網(wǎng)格點,地震記錄總時長9.2s,為了滿足波場延拓穩(wěn)定性條件,選取時間步長0.5ms,需要延拓18401次,保存所有震源波場需要內(nèi)存約246.9GB。因此,SHI等[36]考慮存儲較少的波場信息,進而重建震源波場;王保利等[26]利用存儲若干PML吸收層信息實現(xiàn)了波場重建并將該方法應用到聲波逆時偏移成像中。對于TTI介質(zhì),逆時偏移成

像對數(shù)據(jù)存儲需求更高,為了增強方法的實用性,本文中提出了實用的TTI介質(zhì)的checkpoint波場重建方法。

假定地震波場水平方向為Nx個網(wǎng)格,垂直方向為Nz個網(wǎng)格,為了避免邊界反射噪聲,在有效波場外側(cè)添加PML吸收衰減層(圖1中綠色區(qū)域),在有效波場內(nèi)部設置內(nèi)邊界層(圖1中藍色區(qū)域)。內(nèi)邊界層數(shù)設定為有限差分階數(shù)的一半,如果采用空間12階差分,則內(nèi)邊界為6層。各向異性介質(zhì)checkpoint波場重建由兩次震源波場正傳和一次波場反傳組成。其中一次震源波場正傳在最小時刻和最大時刻之間進行,而另一次震源波場正傳和一次震源波場反傳逐個在兩個checkpoint之間進行。具體如下。

圖1　波場存儲示意

假設震源波場正向傳播的時間點為N×m,N為在此傳播過程中確定的checkpoint個數(shù),m為每兩個checkpoint間的時間點數(shù)。在波場正向延拓過程中,保存每個checkpoint及其前一時間點處的p和q全部波場(圖1中PML吸收層、內(nèi)邊界層和內(nèi)部區(qū))至計算機內(nèi)存,圖2中黑點表示存儲p和q全部波場的時間位置。例如,當波場延拓至checkpointN時,記錄(N×m-1)和N×m時間點對應的p和q全部波場。

以checkpoint(N-1)至checkpointN期間的波場重建為例,需從計算機內(nèi)存中讀取(N-1)×m-1及(N-1)×m時間點對應的p和q全部波場信息作為初值進行m次正向延拓,并將每次延拓波場的內(nèi)邊界層(圖1中藍色區(qū)域)信息保存至GPU的顯存。利用公式(13)將N×m-1和N×m節(jié)點對應的p和q全部波場作為反向延拓的初值,每延拓一次將上一步驟中保存至顯存上p和q的內(nèi)邊界層(圖1中藍色區(qū)域)信息載入到重建波場的對應位置,延拓m次便可重建出checkpoint(N-1)至checkpointN之間對應的震源波場。其它任兩個checkpoint之間波場重建及內(nèi)邊界層的保存辦法同上,進而重建出N×m至最小時刻的震源波場。本文提出的TTI介質(zhì)逆時偏移成像方法是在checkpoint框架下,通過保存內(nèi)邊界層信息實現(xiàn)震源波場重建,而非保存經(jīng)過吸收衰減后的PML吸收層波場信息,因此確保了波場重建精度。另外,在波場重建過程中,與保存所有波場相比,僅保存兩個checkpoint之間波場的內(nèi)邊界層信息,雖然多了一次波場正向延拓的計算量,但是GPU/CPU協(xié)同并行計算技術(shù)在本文方法中的應用,較大地提高了計算效率,充分補償了由此帶來的計算負擔,達到了以計算換存儲的目的。需要指出的是在延拓時間點一定的情況下,checkpoint越多,GPU顯存占用就越少,雖然降低了對GPU硬件的要求,但是增加了計算機內(nèi)存與顯存之間數(shù)據(jù)交換的次數(shù),這將導致計算耗時明顯增長,所以在確保兩checkpoint間的波場內(nèi)邊界占用空間不超過GPU顯存的情況下,可適當減少checkpoint數(shù),以降低數(shù)據(jù)交換頻率,從而達到計算效率的優(yōu)化。

圖2　checkpoint存儲示意

3　各向異性介質(zhì)模型數(shù)據(jù)測試

3.1　TTI介質(zhì)正演數(shù)值模擬

為了驗證上述推導帶有PML吸收邊界qP波方程有限差分格式的正確有效性,測試了均勻TTI介質(zhì)的地震波場傳播。假定模型在水平和垂直方向的網(wǎng)格數(shù)均為320,水平方向和垂直方向的網(wǎng)格間距均為6.25m,沿對稱軸方向的縱波速度vPz=3000m/s,Thomsen各向異性參數(shù)δ=0.1,ε=0.3,PML吸收邊界層數(shù)為50,震源位于模型的中央,主頻為30Hz的雷克子波作為激發(fā)震源,時間采樣間隔為0.5ms。

利用本文推導的公式(12)進行數(shù)值模擬得到的波場快照如圖3所示。圖3a為對稱軸傾角θ=0,沿著對稱軸方向的橫波速度為0,波傳播時間為200ms時p的波場快照。圖3b為對稱軸傾角θ=0,沿著對稱軸方向的橫波速度為0,波傳播時間為350ms時p的波場快照,觀察圖3b可知,PML邊界對地震邊界波場的吸收效果較為理想。圖3c為對稱軸傾角θ=30°,沿著對稱軸方向的橫波速度為0,波傳播時間為200ms時p的波場快照。從圖3c可以看出,qSV波附近位置出現(xiàn)嚴重頻散現(xiàn)象。圖3d為對稱軸傾角θ=30°,沿著對稱軸方向的橫波速度為1500m/s,波傳播時間為200ms時的p的波場快照,可見,加入適當?shù)臋M波速度或者引入σ參數(shù)后,頻散得到了有效抑制,TTI介質(zhì)數(shù)值模擬精度較高。圖3a,圖3b和圖3c 中雖然橫波沿對稱軸方向的速度為0,但是沿著其它方向的橫波速度并不為0,因此出現(xiàn)了qSV波,如圖中菱形狀波場所示。

圖3　p的波場快照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

測試結(jié)果表明,對于圖3d所示的數(shù)值模擬,采用CPU算法共耗時7min,采用CPU/GPU協(xié)同并行加速算法耗時8.5s,計算效率提高約50倍。如果將模型的橫、縱向網(wǎng)格點數(shù)由320增加至3200,數(shù)據(jù)量增加100倍,采用CPU/GPU協(xié)同并行加速算法耗時40.9s,計算耗時是原數(shù)據(jù)的5倍。由此可見,在波場正演模擬中引入CPU/GPU協(xié)同并行加速技術(shù),可以大幅提升計算效率,隨著模型網(wǎng)格點數(shù)的增大,并行計算耗時也相應地增加,但并不是線性增加。

3.2　TTI介質(zhì)逆時偏移試算

利用BP2007復雜構(gòu)造模型檢驗本文介紹的TTI介質(zhì)逆時偏移成像方法。模型水平方向12596個網(wǎng)格點,垂直方向1801網(wǎng)格點,網(wǎng)格間距6.25m,水平長度78.70km,垂直深度11.25km。該模型數(shù)據(jù)共1641炮,每炮800道,最小偏移距37.5m,道間距12.5m,時間方向1151個采樣點,采樣間隔為8ms。本文截取其水平方向43.70km,共626炮數(shù)據(jù)進行TTI介質(zhì)逆時偏移成像測試,截取部分模型參數(shù)如圖4所示。

首先利用截取部分模型來測試基于checkpoint方法重建qP波的可行性和有效性。截取模型的水平方向2000個網(wǎng)格點,垂直方向1801網(wǎng)格點,時間采樣間隔為0.5ms,利用主頻為30Hz的雷克子波作為激發(fā)震源,激發(fā)點位于水平方向1000網(wǎng)格點、垂直方向900網(wǎng)格點處。震源波場傳播至1.5s時的波場快照如圖5a所示,震源波場傳播至最大時刻后,再根據(jù)本文提出的checkpoint方法和內(nèi)邊界保存技術(shù)對震源波場進行分時段的正向傳播和逆向傳播,逆向傳播至1.5s時的震源波場如圖5b所示。正向傳播的震源波場和逆向傳播的震源重建波場在1.5s時刻的波場誤差如圖5c所示,誤差波場振幅與有效波場振幅相差5個數(shù)量級,可忽略不計。由此可知,在checkpoint方法框架下保存波場內(nèi)邊界的方法,可以有效地進行波場重建,豐富了逆時偏移算法的研究,為TTI介質(zhì)逆時偏移方法的工業(yè)化應用帶來了更廣闊的發(fā)展空間。

利用截取的BP2007模型進行偏移試算。單炮偏移區(qū)域水平方向2000個網(wǎng)格點,垂直方向1801個網(wǎng)格點,空間網(wǎng)格間距為6.25m,時間采樣間隔為0.5ms。先將震源子波和地層模型對震源波場正向延拓,保存震源波場部分信息,然后根據(jù)保存的波場信息,基于checkpoint方法對震源波場重建,與此同時將炮集數(shù)據(jù)加載到檢波點得到沿時間逆向的檢波點波場,最后經(jīng)震源歸一化成像條件得到單炮偏移剖面。經(jīng)單炮逆時偏移算法測試,采用CPU/GPU協(xié)同并行加速算法耗時39min,從上述波場重建方法可知該時長受checkpoint數(shù)的影響,而采用CPU算法耗時1755min,GPU硬件的應用可大幅提升逆時偏移的計算效率。對單炮偏移剖面利用拉普拉斯方法去除低頻噪聲后,將626炮偏移結(jié)果進行多炮疊加得到最后的偏移剖面,如圖6所示。從圖6可知,基于qP波方程得到的剖面同相軸連續(xù)性好,兩個陡傾角鹽丘與模型位置匹配。

圖4　TTI介質(zhì)模型參數(shù)a 縱波沿著對稱軸的速度; b Thomsen參數(shù)θ; c Thomsen參數(shù)δ; d Thomsen參數(shù)ε

為了便于分析本文逆時偏移方法的精確性,根據(jù)偏移速度場,結(jié)合常密度參數(shù),計算得到反射系數(shù)模型,如圖7所示。分別在圖6和圖7的水平方向11.25km位置抽取單道數(shù)據(jù)如圖8所示(圖中橫軸為剖面的垂直深度,縱軸是經(jīng)歸一化的單道振幅信息)。由圖8可知,逆時偏移剖面包含了強反射系數(shù)信息,同時弱反射系數(shù)也得以體現(xiàn)。因此,本文的逆時偏移方法可以實現(xiàn)對地下各向異性構(gòu)造體的精細刻畫。

圖5　震源波場重建a 正傳震源波場; b 重建震源波場; c 重建震源波場誤差

圖6　逆時偏移剖面

圖7　反射系數(shù)模型

圖8　水平方向11.25km處單道數(shù)據(jù)對比

4　結(jié)束語

從精準相速度公式出發(fā),推導了基于有限差分數(shù)值解法的TTI介質(zhì)qP波方程的PML吸收邊界格式。由波場快照可知,當波場能量傳播至邊界層后吸收效果良好。在波動方程中引入橫波速度,數(shù)值試驗結(jié)果表明,該方法有效地壓制了qSV波頻散,能夠解決波場傳播不穩(wěn)定的問題。為了解決算法無法在有限顯存GPU端運行的問題,本文在各向異性介質(zhì)偏移成像過程中應用了基于checkpoint方法和保存內(nèi)邊界波場的存儲策略,極大地降低了存儲,雖然該策略增加了一次波場正向延拓計算量以及硬件之間數(shù)據(jù)的傳輸時間,但是偏移成像過程中應用了GPU/CPU并行加速技術(shù),總體上大幅度提升了計算效率,實現(xiàn)以計算換存儲的目的,從而降低對計算硬件的依賴,為TTI介質(zhì)逆時偏移成像的工業(yè)化應用開辟了有效的思路。

[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]秦海旭,吳國忱.TTI介質(zhì)彈性波隨機邊界逆時偏移的實現(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]張千祥,王德利,周進舉.二維TTI介質(zhì)的純P波波動方程數(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ì)有限差分逆時偏移的穩(wěn)定性探討[J].地球物理學報,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]王保利,高靜懷,陳文超,等.地震疊前逆時偏移的有效邊界存儲策略[J].地球物理學報,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]劉國峰,劉洪,李博,等.山地地震資料疊前時間偏移方法及其GPU實現(xiàn)[J].地球物理學報,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]李博,劉紅偉,劉國峰,等.地震疊前逆時偏移算法的CPU/GPU實施對策[J].地球物理學報,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]劉紅偉,李博,劉洪,等.地震疊前逆時偏移高階有限差分算法及GPU實現(xiàn)[J].地球物理學報,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]柯璇,石穎,張瑩瑩,等.地震疊前逆時偏移衰減隨機邊界條件研究[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)格有限差分數(shù)值模擬方法[J].地球物理學報,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