擬空間域彈性波方程交錯網(wǎng)格有限差分逆時偏移

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

1 山東科技大學，海洋科學與工程學院，青島 266590 2 青島海洋科學與技術國家實驗室，海洋地質過程與環(huán)境功能實驗室，青島 266061 3 中國海洋大學，海洋地球科學學院，青島 266100 4 自然資源部，國家深?；毓芾碇行?，青島 266237

0 引言

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

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

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

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

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

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

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

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

1 擬空間域一階速度-應力彈性波方程

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

（1）

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

（2）

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

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

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

（3）

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

（4）

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

（5）

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

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

2.1 擬空間域采樣間隔計算

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

圖3 網(wǎng)格模型中的任意一點P（i,j）周圍有四個擬空間域采樣間隔Fig.3 Four pseudo-space intervals around the point P（i,j）

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

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

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

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

（6）

（7）

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

（8）

（9）

（10）

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

（11）

（12）

（12′）

其中k表示離散的時間點，滿足t=kΔt（Δt表示離散時間步長）.

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

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

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

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

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

（13）

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

（14）

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

（15）

進而可得關于τx的二階導數(shù)表達式：

（16）

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

（18）

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

（19）

此外，關于時間t的二階導數(shù)表達式可寫為：

（20）

（21）

（22）

（23）

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

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

3 模型實驗

3.1 歸一化互相關成像條件

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

（24）

（25）

式中（vp）F為正時縱波波場，（vp）R為逆時縱波波場，（vs）R為逆時橫波波場.

3.2 傾斜界面模型逆時偏移

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

圖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.

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

采用縱波震源激發(fā)，而模型橫波速度由縱波速度和泊松比值計算得到（本實驗所用泊松比為0.25，密度為2000 kg·m-3）.正演模擬時采用主頻為35 Hz的Ricker子波.為保證合成炮集記錄的精度，模擬時縱橫向網(wǎng)格間距均取為1 m，差分精度為時間二階空間十六階，共合成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)格間距模型分別進行常規(guī)彈性波方程和擬空間域彈性波方程有限差分（差分精度為時間二階空間/擬空間十六階）逆時偏移處理.圖8顯示了在0.9 s時刻第76炮的正時波場波前快照，圖9則為獲得的逆時偏移剖面.

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

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

3.3 Marmousi模型逆時偏移

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

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

圖8 逆時偏移正時波場波前快照（0.9 s時刻，第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 逆時偏移剖面（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 逆時偏移地震剖面局部放大圖示（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模型逆時偏移正時波場波前快照（1.9 s時刻，第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模型彈性波逆時偏移剖面（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模型彈性波逆時偏移地震剖面局部放大圖示（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可見，與常規(guī)彈性波方程逆時偏移相比，擬空間域彈性波方程逆時偏移波場中層間反射波明顯減弱（如圖12中紅色橢圓區(qū)域所示），這同樣顯示出擬空間域彈性波方程在復雜構造模型波場延拓時壓制層間反射波的有效性.

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

4 結論與展望

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

（1） 本文推出的高階有限差分格式能夠實現(xiàn)擬空間域彈性波方程的高精度波場延拓，并且其適用于復雜構造模型.

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

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