張曉波，宋鵬，劉保華，譚軍，李金山

1 山東科技大學(xué)，海洋科學(xué)與工程學(xué)院，青島 266590 2 青島海洋科學(xué)與技術(shù)國(guó)家實(shí)驗(yàn)室，海洋地質(zhì)過(guò)程與環(huán)境功能實(shí)驗(yàn)室，青島 266061 3 中國(guó)海洋大學(xué)，海洋地球科學(xué)學(xué)院，青島 266100 4 自然資源部，國(guó)家深海基地管理中心，青島 266237

0 引言

地震勘探是目前探測(cè)復(fù)雜地下地質(zhì)構(gòu)造最重要的地球物理方法之一，主要包括地震信息采集、數(shù)據(jù)處理和成像解釋三大環(huán)節(jié)（Yilmaz,2001）.而地震數(shù)據(jù)處理是連接信息采集與成像解釋環(huán)節(jié)的重要橋梁.長(zhǎng)期以來(lái)，精確的地下構(gòu)造偏移成像一直是地震數(shù)據(jù)處理中的研究熱點(diǎn)（Guo et al.,2002；Chattopadhyay and McMechan,2008）.

基于雙程波動(dòng)方程理論進(jìn)行逆時(shí)偏移成像是當(dāng)前公認(rèn)的精確構(gòu)造成像方法（Mora,1987;Sun and McMechan,2001;Sirgue et al.,2010;Plessix and Perkins,2010;王保利等,2012;Moradpouri et al.,2017）.逆時(shí)偏移（Reverse Time Migration,RTM）算法出現(xiàn)于20世紀(jì)80年代早期（Whitmore,1983;McMechan,1983;Chang and McMechan,1987,1990,1994），其不需要對(duì)波動(dòng)方程進(jìn)行分解，避免了因?qū)Σ▌?dòng)方程近似而造成的傾角限制，因此該方法能夠適用于復(fù)雜構(gòu)造成像，一直以來(lái)都是偏移成像領(lǐng)域的重要研究?jī)?nèi)容（Clapp,2009;劉紅偉等,2010;劉欣欣等,2013;石穎等,2015;宋鵬等,2015;李慶洋等,2017;Li et al.,2018,2020）.

近年來(lái)，國(guó)內(nèi)外地震勘探技術(shù)發(fā)展迅速，特別是多波多分量地震數(shù)據(jù)采集技術(shù)已日漸成熟并且在實(shí)際油氣勘探中取得了良好的應(yīng)用效果，值得大力推廣（張樹(shù)林和張懿,2014）.為了充分利用多波多分量地震數(shù)據(jù)中豐富的縱、橫波信息，實(shí)現(xiàn)地下構(gòu)造精確成像，大力發(fā)展彈性波方程逆時(shí)偏移已勢(shì)在必行.目前一些學(xué)者已開(kāi)展相關(guān)研究工作，如：Chang和McMechan首先將二維逆時(shí)偏移推廣到彈性波領(lǐng)域（Chang and McMechan,1987），隨后又將其推廣至三維情況（Chang and McMechan,1990,1994）；Luo 等（2013）利用譜元法和基于伴隨矩陣的成像條件實(shí)現(xiàn)了彈性波逆時(shí)偏移；王玉鳳（2014）基于SEG/EAGE模型實(shí)現(xiàn)了多波多分量數(shù)據(jù)彈性波方程逆時(shí)偏移；Duan和Sava（2017）實(shí)現(xiàn)了彈性波方程逆時(shí)偏移成像道集的角度域分解；Du等（2014）和Gong等（2018）分別基于不同的縱橫波分量實(shí)現(xiàn)了彈性波逆時(shí)偏移成像；Zhang和Shi（2019）研究了彈性波逆時(shí)偏移的成像條件；Qu 等（2020）為解決常規(guī)縱橫波波場(chǎng)分離方法存在的PS和SP分量動(dòng)力學(xué)特性改變以及極性反轉(zhuǎn)等問(wèn)題，提出了一種基于曲線坐標(biāo)的OBC數(shù)據(jù)波場(chǎng)分離彈性逆時(shí)偏移方法.

隨著研究的不斷深入，彈性波方程逆時(shí)偏移方法已日趨完善，然而在當(dāng)前的實(shí)現(xiàn)和應(yīng)用過(guò)程中，其仍存在以下兩個(gè)方面的問(wèn)題有待進(jìn)一步研究.首先，波動(dòng)方程有限差分逆時(shí)偏移技術(shù)在當(dāng)前的實(shí)現(xiàn)中均普遍應(yīng)用到重要的數(shù)值計(jì)算方法——有限差分法（Bartolo et al.,2017;Ren and Li,2017；趙明哲等,2022）.眾所周知，傳統(tǒng)的波動(dòng)方程有限差分計(jì)算方法本身存在著難以克服的固有問(wèn)題——空間域的矩形網(wǎng)格剖分必然造成地震速度界面的畸變（劉立彬等,2020）.圖1給出了一個(gè)包含傾斜速度界面模型的常規(guī)有限差分法網(wǎng)格剖分示意圖.圖中黑色實(shí)線為傾斜速度界面，而紅色實(shí)心點(diǎn)和藍(lán)色實(shí)心點(diǎn)分別表示界面兩側(cè)的速度.顯然，矩形網(wǎng)格剖分的結(jié)果使得原本光滑的傾斜速度界面畸變?yōu)槊黠@的階梯狀折線，當(dāng)?shù)卣鸩▊鞑ブ吝@樣的界面時(shí)會(huì)形成一系列的假散射（褚春雷和王修田,2005）；此外地震波在遇到速度分界面時(shí)會(huì)產(chǎn)生大量的層間反射波，其在偏移剖面上形成較強(qiáng)的低波數(shù)干擾，其會(huì)嚴(yán)重影響剖面的信噪比并可產(chǎn)生偏移假象（Liu et al.,2011a），從而影響偏移剖面的精度.這兩個(gè)問(wèn)題的有效解決對(duì)提高逆時(shí)偏移的成像質(zhì)量具有重要意義.

圖1 常規(guī)有限差分法網(wǎng)格剖分示意圖Fig.1 Schematic diagram of regular meshing in conventional finite difference method

為了消除逆時(shí)偏移過(guò)程中的假散射和界面畸變的現(xiàn)象，一些學(xué)者采用空間可變網(wǎng)格（朱生旺和魏修成,2005；朱生旺等,2007;黃超和董良國(guó),2009；姜占東等,2021）的思路，在介質(zhì)變化劇烈的區(qū)域采用精細(xì)網(wǎng)格，在介質(zhì)變化平緩的區(qū)域采用較粗網(wǎng)格剖分，但該方法仍沒(méi)有擺脫矩形網(wǎng)格的局限性；褚春雷和王修田（2005）根據(jù)有限元的思想采用基于非規(guī)則三角形網(wǎng)格的有限差分方法實(shí)現(xiàn)有限差分模擬，該方法可以更精細(xì)地描述起伏界面，但是其計(jì)算量與傳統(tǒng)矩形網(wǎng)格有限差分相比明顯增加.

而對(duì)于由于層間反射波引起的低波數(shù)噪聲干擾壓制問(wèn)題，當(dāng)前主要有如下幾類方法：第一類是背向反射壓制方法，其通常是基于能夠弱化反射波場(chǎng)的波動(dòng)方程進(jìn)行波場(chǎng)延拓，以達(dá)到壓制層間反射波、減弱低波數(shù)干擾的目的.Baysal等（1984）首先提出了基于常波阻抗假設(shè)的無(wú)反射聲波方程，可顯著壓制垂向附近入射的地震波背向反射；宋鵬（2005）和Zhang 等（2010）改進(jìn)了無(wú)反射聲波方程，提升了背向反射的壓制效果.第二類方法為濾波類方法，如高通濾波、拉普拉斯濾波等.Mulder和Plessix（2004）直接采用高通濾波方法對(duì)成像剖面進(jìn)行去噪處理；Zhang和Sun（2009）采用Laplacian濾波方法對(duì)常規(guī)逆時(shí)偏移結(jié)果進(jìn)行濾波.第三類是行波分離方法，該方法對(duì)震源波場(chǎng)和檢波點(diǎn)波場(chǎng)分離成不同方向的行波，然后提取有效波場(chǎng)分量參與成像，從而實(shí)現(xiàn)地下構(gòu)造的精確成像.Liu 等（2011a）首先基于F-K變換實(shí)現(xiàn)了波場(chǎng)分離逆時(shí)偏移成像，有效地壓制了偏移噪聲；Fei等（2015）、王一博等（2016）應(yīng)用Hilbert變換實(shí)現(xiàn)了行波分離的逆時(shí)偏移成像；Chen和He（2014）采用Poynting矢量實(shí)現(xiàn)了震源波場(chǎng)與檢波點(diǎn)波場(chǎng)上、下、左、右四個(gè)方向的行波分離，有效提高低波數(shù)噪聲的壓制效果.以上三類方法各有優(yōu)劣，相對(duì)而言，背向反射壓制與波場(chǎng)分離這兩類方法是在成像過(guò)程中壓制低波數(shù)噪聲，其對(duì)于最終成像質(zhì)量的提高更為有效，對(duì)于這二者來(lái)講，背向反射壓制方法更具計(jì)算效率優(yōu)勢(shì)；而濾波類方法雖可有效壓制低波數(shù)噪聲，但該類方法是在成像計(jì)算之后壓制低波數(shù)噪聲，對(duì)于復(fù)雜構(gòu)造模型其成像剖面的精度往往受到一定影響，且該類方法還存在濾波器閾值范圍選擇困難，易損傷有效信息，有時(shí)還會(huì)引入大量高頻噪聲（陳康和吳國(guó)忱,2012）等問(wèn)題，因此該類方法在實(shí)際數(shù)據(jù)處理時(shí)需謹(jǐn)慎使用，或與其他類方法搭配使用以期達(dá)到更優(yōu)的成像效果.

2005年SEG年會(huì)上，Wang 等（2005）通過(guò)將傳統(tǒng)的波動(dòng)方程從時(shí)間-空間域變換到時(shí)間-走時(shí)域（或稱作“走時(shí)域”），推導(dǎo)出了一種擬空間域聲波方程，其既可有效解決常規(guī)波動(dòng)方程逆時(shí)偏移中的彎曲界面假散射和界面畸變的問(wèn)題，同時(shí)像無(wú)反射聲波方程類似，該方程還可顯著減少波場(chǎng)延拓中的背向反射，達(dá)到壓制成像剖面中的低波數(shù)噪聲，提高成像質(zhì)量的目的.但Wang等（2005）僅實(shí)現(xiàn)了基于擬空間域二階差分精度有限差分模擬，其模擬精度較低，限制了該方法在實(shí)際中的應(yīng)用.本文在詳細(xì)討論擬空間域彈性波方程原理的基礎(chǔ)上，推導(dǎo)了任意偶數(shù)階擬空間域一階速度-應(yīng)力彈性波方程交錯(cuò)網(wǎng)格有限差分格式，并給出了其穩(wěn)定性條件，實(shí)現(xiàn)了高精度的擬空間域一階速度-應(yīng)力彈性波方程交錯(cuò)網(wǎng)格有限差分逆時(shí)偏移.

1 擬空間域一階速度-應(yīng)力彈性波方程

在二維各向同性介質(zhì)中，拉梅常數(shù)λ和μ與介質(zhì)縱波速度vp、橫波速度vs以及密度ρ之間滿足關(guān)系：

（1）

根據(jù)式（1）和常規(guī)一階速度-應(yīng)力彈性波方程（Virieux,1986;Coutant et al.,1995;Liu and Sen,2011b），二維各向同性彈性波方程可表示成：

（2）

其中，x、z分別代表空間的橫向坐標(biāo)和縱向坐標(biāo)，σp表示應(yīng)力波場(chǎng)的縱波分量，σxxs表示水平方向上正應(yīng)力的橫波分量，σzzs表示垂直方向上正應(yīng)力的橫波分量，σxz表示切應(yīng)力波場(chǎng)，vx表示水平速度分量波場(chǎng)，vz表示垂直速度分量波場(chǎng)，s（t）為震源函數(shù)，t表示時(shí)間.

圖2 空間模型離散化之后的差分網(wǎng)格劃分示意圖Fig.2 Illustration of the spatial model after discretization

當(dāng)基于式（2）所示的彈性波方程進(jìn)行有限差分波場(chǎng)延拓時(shí)，將時(shí)間和空間離散化后（如圖2所示），令i和j分別代表x和z方向上網(wǎng)格點(diǎn)坐標(biāo).假設(shè)空間無(wú)限小的網(wǎng)格長(zhǎng)度為Δξ（ξ可以代表x或z），地震波在該網(wǎng)格上縱橫波走時(shí)分別為Δτ（p）ξ和Δτ（s）ξ，則空間網(wǎng)格Δξ與走時(shí)Δτ（p）ξ和Δτ（s）ξ之間滿足關(guān)系式：

（3）

根據(jù)式（3）可以將應(yīng)力σ（σ代表σp,σxxs,σzzs或σxz）以及速度分量vx、vz對(duì)空間的導(dǎo)數(shù)變換為：

（4）

然后將式（4）代入式（2）中可得：

（5）

式（5）即為擬空間域一階速度-應(yīng)力彈性波方程.

2 擬空間域一階速度-應(yīng)力彈性波方程數(shù)值計(jì)算

2.1 擬空間域采樣間隔計(jì)算

通常情況下，為了采用有限差分方法描述擬空間域彈性波方程（如式（5）所示），首先要按照傳統(tǒng)的有限差分方法將連續(xù)的空間模型變成離散化的空間網(wǎng)格模型，然后在每個(gè)空間網(wǎng)格線Δξ上，根據(jù)其對(duì)應(yīng)的網(wǎng)格速度值計(jì)算縱橫波走時(shí)Δτ（p）ξ和Δτ（s）ξ.為了便于描述，本節(jié)中將Δτ（p）ξ和Δτ（s）ξ統(tǒng)稱為擬空間域采樣間隔Δτ（φ）ξ，其中φ表示縱波p或橫波s.

圖3 網(wǎng)格模型中的任意一點(diǎn)P（i,j）周?chē)兴膫€(gè)擬空間域采樣間隔Fig.3 Four pseudo-space intervals around the point P（i,j）

圖4 傾斜速度界面的擬空間域采樣間隔計(jì)算示意圖Fig.4 Partial schematic illustration of a mesh model after the regular meshing of a velocity interface model including a slanted interface

理論上講，在擬空間域不再存在速度界面的畸變問(wèn)題，甚至還會(huì)減弱相鄰網(wǎng)格點(diǎn)之間模型參量的突變程度，從而可在偏移計(jì)算中期望降低假散射和界面反射.

2.2 擬空間域彈性波方程2N階精度有限差分格式

基于擬空間域彈性波方程進(jìn)行有限差分?jǐn)?shù)值模擬時(shí)，為了提高模擬的精度，減少數(shù)值頻散的影響，需要提高差分的精度，因此本文推導(dǎo)了擬空間域彈性波方程2N階精度交錯(cuò)網(wǎng)格有限差分表達(dá)式.

（6）

（7）

1,2,…,N-1,N）），在τ（p）x=τi+1/2處的2N階Taylor級(jí)數(shù)展開(kāi)式為：

（8）

（9）

（10）

將式（10）求得的差分系數(shù)代入式（9）可以得到σp對(duì)變量τ（p）x在（τi+1/2,τj）處的一階導(dǎo)數(shù)2N階精度差分表達(dá)式：

（11）

（12）

（12′）

其中k表示離散的時(shí)間點(diǎn)，滿足t=kΔt（Δt表示離散時(shí)間步長(zhǎng)）.

在中心波場(chǎng)計(jì)算區(qū)域，應(yīng)用（12）式即可實(shí)現(xiàn)一階速度-應(yīng)力彈性波方程的時(shí)間二階擬空間域2N階精度有限差分?jǐn)?shù)值模擬；而在人工邊界區(qū)域，為有效壓制人工邊界反射，還需進(jìn)行吸收邊界處理，本文的數(shù)值實(shí)驗(yàn)均采用完全匹配層（Perfectly Matched Layer，PML）吸收邊界條件（Collino and Tsogka,2001;任志明和劉洋,2014;張曉波等,2016;張曉波,2017）.

2.3 擬空間域一階速度-應(yīng)力彈性波方程穩(wěn)定性條件

結(jié)合常規(guī)一階速度-應(yīng)力彈性波方程交錯(cuò)網(wǎng)格有限差分格式的穩(wěn)定性條件（裴正林和牟永光,2003）以及平面諧波分析的方法，可以推導(dǎo)出一階速度-應(yīng)力彈性波方程擬空間域交錯(cuò)網(wǎng)格有限差分格式的穩(wěn)定性條件.

首先定義擬空間域平面諧波變量u：

u=u0eiω nΔteikτxjΔτxeikτzkΔτz，

（13）

其中u0表示初始波場(chǎng)，ω代表圓頻率，kτx和kτz分別代表關(guān)于τx和τz的波數(shù)，n、j和k分別代表t、τx和τz方向上離散網(wǎng)格點(diǎn)的坐標(biāo)，Δτx和Δτz分別代表τx和τz方向上的擬空間域采樣間隔，e代表自然對(duì)數(shù)的底，i代表虛數(shù)單位.根據(jù)式（13）可以得到：

（14）

將式（14）代入一階導(dǎo)數(shù)差分格式表達(dá)式中為：

（15）

進(jìn)而可得關(guān)于τx的二階導(dǎo)數(shù)表達(dá)式：

（16）

將式（17）中差分系數(shù)分別用關(guān)于τx的縱波和橫波擬空間域差分系數(shù)代替，可得關(guān)于τx的縱波和橫波擬空間域二階導(dǎo)數(shù)表達(dá)式：

（18）

類似地，可推導(dǎo)出關(guān)于τz的縱波和橫波擬空間域二階導(dǎo)數(shù)表達(dá)式：

（19）

此外，關(guān)于時(shí)間t的二階導(dǎo)數(shù)表達(dá)式可寫(xiě)為：

（20）

（21）

（22）

（23）

式中模型最大速度vmax=vp且Δh=Δx=Δz，將式（22）代入式（23）可以得到庫(kù)朗數(shù)隨差分階數(shù)的變化關(guān)系，如圖5所示.從圖中可以看出，庫(kù)朗數(shù)隨著有限差分階數(shù)的增加而減小，并且當(dāng)階數(shù)為2時(shí)，庫(kù)朗數(shù)值小于1，符合數(shù)值計(jì)算的穩(wěn)定性要求.

圖5 庫(kù)朗數(shù)隨有限差分階數(shù)的變化關(guān)系圖Fig.5 Diagram of Courant number variation with finite difference orders

3 模型實(shí)驗(yàn)

3.1 歸一化互相關(guān)成像條件

本文在逆時(shí)偏移中采用歸一化互相關(guān)成像條件（Kaelin and Guitton,2006）實(shí)現(xiàn)偏移成像，其實(shí)現(xiàn)過(guò)程是利用正時(shí)波場(chǎng)的零延遲互相關(guān)的結(jié)果對(duì)正時(shí)波場(chǎng)與逆時(shí)波場(chǎng)零延遲互相關(guān)成像進(jìn)行歸一化處理（如式（24）和式（25）所示）.這里根據(jù)縱橫波波場(chǎng)分離方法（李振春等,2007），對(duì)式（5）進(jìn)行分解可以得到縱波波場(chǎng)vp={vxp,vzp}和橫波波場(chǎng)vs={vxs,vzs}，其中vxp和vxs分別為水平速度分量vx的縱波波場(chǎng)和橫波波場(chǎng)；vzp和vzs分別為垂直速度分量vz的縱波波場(chǎng)和橫波波場(chǎng).然后將縱橫波分離后的波場(chǎng)代入式（24）和式（25）中，即可求出正時(shí)縱波波場(chǎng)和逆時(shí)縱波波場(chǎng)（PP 波）成像結(jié)果IPP和正時(shí)縱波波場(chǎng)和逆時(shí)橫波波場(chǎng)（PS波）成像結(jié)果IPS.式（24）和式（25）為：

（24）

（25）

式中（vp）F為正時(shí)縱波波場(chǎng)，（vp）R為逆時(shí)縱波波場(chǎng)，（vs）R為逆時(shí)橫波波場(chǎng).

3.2 傾斜界面模型逆時(shí)偏移

本實(shí)驗(yàn)的主要目的是檢驗(yàn)擬空間域聲波方程逆時(shí)偏移在解決速度界面畸變以及壓制界面假散射和層間反射波方面的有效性.

圖6 含傾斜界面的兩層速度模型（a） 原光滑界面模型；（b） 10 m網(wǎng)格間距的模型.Fig.6 Two-layer velocity model with a slanted interface（a） Original model with a smooth slanted interface;（b） Model with 10 m grid interval.

實(shí)驗(yàn)采用含有傾斜界面的兩層速度模型，如圖6a所示，其橫向和縱向長(zhǎng)度分別為4000 m和2000 m，界面上下兩側(cè)縱波速度分別為2500 m·s-1和3500 m·s-1，密度為2000 kg·m-3.將此傾斜界面模型以10 m的縱橫向網(wǎng)格間距剖分后所得網(wǎng)格模型如圖6b所示，可見(jiàn)原光滑的速度界面已變?yōu)槊黠@的階梯狀界面（見(jiàn)圖中白色箭頭所指位置）.實(shí)驗(yàn)中建立道固定、炮移動(dòng)的觀測(cè)系統(tǒng)，炮點(diǎn)位于500 m至3480 m之間，炮間隔為20 m，共150炮；每炮接收道數(shù)為401道，各接收道位于0 m至4000 m之間，道間隔為10 m；炮點(diǎn)與接收點(diǎn)深度均為10 m.

采用縱波震源激發(fā)，而模型橫波速度由縱波速度和泊松比值計(jì)算得到（本實(shí)驗(yàn)所用泊松比為0.25，密度為2000 kg·m-3）.正演模擬時(shí)采用主頻為35 Hz的Ricker子波.為保證合成炮集記錄的精度，模擬時(shí)縱橫向網(wǎng)格間距均取為1 m，差分精度為時(shí)間二階空間十六階，共合成150炮地震記錄，其中第76炮記錄如圖7所示.

圖7 合成炮集記錄（第76炮）（a） 水平速度分量；（b） 垂直速度分量.Fig.7 Synthetic shot gather record example （76th shot gather）（a） Horizontal velocity component；（b） Vertical velocity component.

基于10 m網(wǎng)格間距模型分別進(jìn)行常規(guī)彈性波方程和擬空間域彈性波方程有限差分（差分精度為時(shí)間二階空間/擬空間十六階）逆時(shí)偏移處理.圖8顯示了在0.9 s時(shí)刻第76炮的正時(shí)波場(chǎng)波前快照，圖9則為獲得的逆時(shí)偏移剖面.

從圖8可以看出，常規(guī)彈性波方程的波場(chǎng)中存在明顯的界面假散射（如圖8a、c中橢圓區(qū)域所示）而擬空間域彈性波方程的波場(chǎng)中并無(wú)明顯的界面假散射（如圖8b、d中橢圓內(nèi)區(qū)域所示）.對(duì)比圖8中箭頭處界面反射波可以看出，擬空間域彈性波方程對(duì)界面反射波（尤其是近垂直入射的反射波）壓制明顯.由此可以證明擬空間域彈性波方程在壓制界面假散射和層間反射波方面的有效性.

為更直觀的對(duì)比兩種方法偏移剖面中傾斜界面的形態(tài)，對(duì)圖9中橢圓區(qū)域內(nèi)界面同相軸進(jìn)行放大顯示，如圖10所示.可以看出，常規(guī)彈性波方程逆時(shí)偏移剖面中傾斜界面形態(tài)（圖10a、c中紅色虛線）相比于真實(shí)界面形態(tài)（圖10a、c紅色實(shí)線）出現(xiàn)明顯畸變，而擬空間域彈性波方程逆時(shí)偏移剖面中傾斜界面的形態(tài)則與真實(shí)界面形態(tài)（圖10b、d中紅色實(shí)線）基本吻合.由此說(shuō)明了擬空間域彈性波方程逆時(shí)偏移在解決界面畸變問(wèn)題方面的有效性.

3.3 Marmousi模型逆時(shí)偏移

Marmousi模型是包含大量速度界面、陡傾角構(gòu)造以及劇烈速度變化的復(fù)雜構(gòu)造網(wǎng)格速度模型，模型橫向和縱向長(zhǎng)度分別為9200 m和3000 m，橫縱向網(wǎng)格間距分別為5 m和4 m，模型橫波速度由縱波速度和泊松比值計(jì)算得到（本實(shí)驗(yàn)所用泊松比為0.25，密度為2000 kg·m-3），縱波和橫波速度模型如圖11a、b所示.實(shí)驗(yàn)中采用右邊放炮、左邊接收的單邊觀測(cè)系統(tǒng)，共426炮，每炮104道接收，炮間隔和道間隔均為25 m，炮點(diǎn)和接收道深度均為8 m.采用縱波震源激發(fā)，震源子波采用主頻為35 Hz的Ricker子波.通過(guò)常規(guī)彈性波方程有限差分方法（差分精度為時(shí)間二階空間八階）正演模擬合成426炮地震記錄.

基于合成炮集記錄，分別進(jìn)行常規(guī)彈性波方程和擬空間域彈性波方程有限差分（差分精度為時(shí)間二階空間/擬空間八階）逆時(shí)偏移處理.圖12顯示了1.9 s時(shí)刻、第138炮的正時(shí)波場(chǎng)波前快照，圖13則給出了相應(yīng)的彈性波逆時(shí)偏移剖面.

圖8 逆時(shí)偏移正時(shí)波場(chǎng)波前快照（0.9 s時(shí)刻，第76炮）（a） 基于常規(guī)彈性波方程的垂直速度分量;（b） 基于擬空間域彈性波方程的垂直速度分量;（c） 基于常規(guī)彈性波方程的水平速度分量;（d） 基于擬空間域彈性波方程的水平速度分量.Fig.8 Snapshot of forward time wavefield in reverse time migration （76th shot at 0.9 s）（a） Vertical velocity component based on the conventional elastic wave equation;（b） Vertical velocity component based on the pseudo-space elastic wave equation;（c） Horizontal velocity component based on the conventional elastic wave equation;（d） Horizontal velocity component based on the pseudo-space elastic wave equation.

圖9 逆時(shí)偏移剖面（a） 基于常規(guī)彈性波方程的IPP剖面；（b） 基于擬空間域彈性波方程的IPP剖面；（c） 基于常規(guī)彈性波方程的IPS剖面；（d） 基于擬空間域彈性波方程的IPS剖面.Fig.9 Profile of reverse time migration（a） IPP based on the conventional elastic wave equation；（b） IPP based on the pseudo-space elastic wave equation；（c） IPS based on the conventional elastic wave equation；（d） IPS based on the pseudo-space elastic wave equation.

圖10 逆時(shí)偏移地震剖面局部放大圖示（a） 基于常規(guī)彈性波方程的IPP剖面；（b） 基于擬空間域彈性波方程的IPP剖面；（c） 基于常規(guī)彈性波方程的IPS剖面；（d） 基于擬空間域彈性波方程的IPS剖面.Fig.10 Local magnification of a reverse time migration profile（a） IPP based on the conventional elastic wave equation；（b） IPP based on the pseudo-space elastic wave equation；（c） IPS based on the conventional elastic wave equation；（d） IPS based on the pseudo-space elastic wave equation.

圖11 Marmousi網(wǎng)格速度模型（a） 縱波速度模型；（b） 橫波速度模型.Fig.11 Grid velocity model of Marmousi（a） Primary velocity model；（b） Shear velocity model.

圖12 Marmousi模型逆時(shí)偏移正時(shí)波場(chǎng)波前快照（1.9 s時(shí)刻，第138炮）（a） 基于常規(guī)彈性波方程的垂直速度分量;（b） 基于擬空間域彈性波方程的垂直速度分量;（c） 基于常規(guī)彈性波方程的水平速度分量;（d） 基于擬空間域彈性波方程的水平速度分量.Fig.12 Forward time wavefield wavefront snapshot in reverse time migration for the Marmousi model （138th shot at 1.9 s in time）（a） Vertical velocity component based on the conventional elastic wave equation;（b） Vertical velocity component based on the pseudo-space elastic wave equation;（c） Horizontal velocity component based on the conventional elastic wave equation;（d） Horizontal velocity component based on the pseudo-space elastic wave equation.

圖13 Marmousi模型彈性波逆時(shí)偏移剖面（a） 基于常規(guī)彈性波方程的IPP剖面；（b） 基于擬空間域彈性波方程的IPP剖面；（c） 基于常規(guī)彈性波方程的IPS剖面；（d） 基于擬空間域彈性波方程的IPS剖面.Fig.13 Reverse time migration profiles for the Marmousi model（a） IPP based on the conventional elastic wave equation；（b） IPP based on the pseudo-space elastic wave equation；（c） IPS based on the conventional elastic wave equation；（d） IPS based on the pseudo-space elastic wave equation.

圖14 Marmousi模型彈性波逆時(shí)偏移地震剖面局部放大圖示（a） 網(wǎng)格模型；（b） 基于常規(guī)彈性波方程；（c） 基于擬空間域彈性波方程.Fig.14 Local magnification of reverse time migration profiles for the Marmousi model （a） Grid model；（b） Based on the conventional elastic wave equation；（c） Based on the pseudo-space elastic wave equation.

由圖12可見(jiàn)，與常規(guī)彈性波方程逆時(shí)偏移相比，擬空間域彈性波方程逆時(shí)偏移波場(chǎng)中層間反射波明顯減弱（如圖12中紅色橢圓區(qū)域所示），這同樣顯示出擬空間域彈性波方程在復(fù)雜構(gòu)造模型波場(chǎng)延拓時(shí)壓制層間反射波的有效性.

為更好地對(duì)比逆時(shí)偏移的效果，將圖13中紅色矩形區(qū)域內(nèi)的局部偏移剖面進(jìn)行放大顯示（如圖14所示）.可以看出，與常規(guī)彈性波方程逆時(shí)偏移相比，擬空間域彈性波方程逆時(shí)偏移所得剖面中構(gòu)造清晰、同相軸連續(xù)性更好（紅色箭頭處所示），從而說(shuō)明擬空間域彈性波方程逆時(shí)偏移的成像質(zhì)量要優(yōu)于常規(guī)彈性波方程逆時(shí)偏移.

4 結(jié)論與展望

本文深入討論了擬空間域彈性波方程的原理，在此基礎(chǔ)上，推導(dǎo)了擬空間域彈性波方程高階有限差分格式，并給出了差分格式的穩(wěn)定性條件，實(shí)現(xiàn)了高精度的擬空間域彈性波方程逆時(shí)偏移.理論分析和模型實(shí)驗(yàn)得出結(jié)論如下：

（1） 本文推出的高階有限差分格式能夠?qū)崿F(xiàn)擬空間域彈性波方程的高精度波場(chǎng)延拓，并且其適用于復(fù)雜構(gòu)造模型.

（2） 若在計(jì)算擬空間采樣間隔時(shí)引入速度界面信息，則擬空間域彈性波方程高階有限差分逆時(shí)偏移能夠避免常規(guī)彈性波方程逆時(shí)偏移中彎曲界面形態(tài)畸變問(wèn)題；此外基于該方法進(jìn)行波場(chǎng)延拓時(shí)可有效壓制彎曲界面的假散射問(wèn)題，并能顯著降低層間反射波，因此可以減少剖面上的偏移假象，從而顯著提高成像的質(zhì)量.

當(dāng)然，目前擬空間域彈性波方程高階有限差分逆時(shí)偏移方法對(duì)于非垂直入射的層間反射波壓制效果尚不理想，因此進(jìn)一步改進(jìn)該方法的層間反射波壓制效果，提高逆時(shí)偏移成像的精度，并將其發(fā)展到三維波動(dòng)方程的逆時(shí)偏移中將是下一步的研究工作.