有限差分的時(shí)間和空間頻散及其對(duì)逆時(shí)偏移和全波形反演的影響

任志明， 戴雪， 包乾宗， 蔡曉慧， 劉洋

1 長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院， 西安 710064 2 南京工業(yè)大學(xué)巖土工程研究所， 南京 210009 3 西安郵電大學(xué)自動(dòng)化學(xué)院， 西安 710121

0 引言

有限差分法計(jì)算量小、容易實(shí)現(xiàn)，是地球物理學(xué)和地震學(xué)領(lǐng)域最常用的波動(dòng)方程數(shù)值解法，也被廣泛應(yīng)用于逆時(shí)偏移和全波形反演中進(jìn)行波場(chǎng)正反向延拓.但由于時(shí)間和空間離散，有限差分法會(huì)產(chǎn)生數(shù)值頻散，不可避免地降低正演、成像和反演的精度.為改善其精度，不同學(xué)者先后提出規(guī)則網(wǎng)格和交錯(cuò)網(wǎng)格(Virieux，1986；Luo and Schuster，1990；董良國(guó)等，2000)、顯式和隱式(Liu and Sen，2009a；Kosloff et al.，2010；Chu and Stoffa，2012)、空間域和時(shí)空域(Liu and Sen，2019b；Ren and Liu，2015；蔡曉慧等，2015)、非規(guī)則網(wǎng)格和無網(wǎng)格(裴正林，2004；孫輝和張劍鋒，2019；Li et al.，2017；Liu et al.，2019；劉立彬等，2020；姜占東等，2021)有限差分.有限差分法在勘探地球物理和計(jì)算地震學(xué)的詳細(xì)進(jìn)展請(qǐng)參考Carcione等(2002)、劉洋(2014)和Moczo 等(2014)，這里不再贅述.

有限差分法的模擬精度由時(shí)間和空間導(dǎo)數(shù)的離散方式?jīng)Q定.提高空間精度的方法主要包括：增大空間差分算子的長(zhǎng)度；優(yōu)化空間差分系數(shù)(Zhang and Yao，2013；Liu，2014；Yang et al.，2017)；偽譜法(Kosloff and Baysal，1982；Fornberg，1987；Liu and Wei，2005；劉財(cái)?shù)龋?013；李青陽(yáng)等，2018).采用長(zhǎng)的差分算子逼近空間偏導(dǎo)數(shù)會(huì)增加計(jì)算量，且存在飽和效應(yīng)(即算子長(zhǎng)度繼續(xù)增大對(duì)模擬精度的改善非常小).常規(guī)有限差分法的差分系數(shù)是基于泰勒級(jí)數(shù)展開得到的，只能保證低波數(shù)區(qū)域的模擬精度.優(yōu)化有限差分法采用優(yōu)化算法(如最小二乘、模擬退火等)逼近頻散關(guān)系求取差分系數(shù)，可以大幅提高中-高波數(shù)區(qū)域的模擬精度.偽譜法通過正反傅里葉變換求取空間偏導(dǎo)數(shù)，理論上具有無窮階/譜精度，但需要更多的計(jì)算時(shí)間.提高有限差分時(shí)間精度的方法包括：Lax-Wendroff法(Dablain，1986；Chen，2009；Amundsen and Pedersen，2017)；新差分模板法(Liu and Sen，2013；Wang et al.，2016；Tang and Huang，2014；張保慶等，2016；Ren et al.，2017)；時(shí)間頻散校正法(Stork，2013；Wang and Xu，2015；Koene et al.，2018).Lax-Wendroff法通過空間偏導(dǎo)數(shù)替換高階時(shí)間偏導(dǎo)數(shù)來改善時(shí)間精度.新差分模板法采用非“十字”或“正交”模板逼近空間偏導(dǎo)數(shù).Liu和Sen(2013)給出一種“菱形”差分模板，新方法可以獲得任意偶數(shù)階時(shí)間精度，但計(jì)算量較大.為平衡精度與效率，Wang等(2016)提出基于“十字+菱形”模板的有限差分法.Tang和Huang(2014)和Ren等(2017)發(fā)展了時(shí)間4～6階和任意偶數(shù)階交錯(cuò)網(wǎng)格有限差分法.新差分模板法的差分系數(shù)隨速度變化，相對(duì)比較復(fù)雜.時(shí)間頻散校正法通過正反時(shí)間頻散變換向觀測(cè)記錄中增加或消除時(shí)間頻散.Wang和Xu(2015)介紹了有限差分的正反時(shí)間頻散變換對(duì)，并將其應(yīng)用于逆時(shí)偏移中消除時(shí)間頻散對(duì)成像的影響.Koene等(2018)指出采用共軛算子代替逆算子的不足，推導(dǎo)了更精確的正時(shí)間頻散算子，使得輸入信號(hào)經(jīng)過正反變換后能夠被精確重建.與其他方法相比，時(shí)間頻散校正法的效率更高，且可以徹底消除波場(chǎng)中的時(shí)間頻散誤差.

波動(dòng)方程正演是逆時(shí)偏移和全波形反演的基本單元，成像和反演的精度依賴于所采用的數(shù)值模擬算法.到目前為止，單獨(dú)研究有限差分的時(shí)間和空間頻散特性及其對(duì)逆時(shí)偏移和全波形反演的影響的文獻(xiàn)相對(duì)較少.本文采用時(shí)間有限差分+空間偽譜法和時(shí)間頻散校正+空間有限差分分離時(shí)間頻散和空間頻散；詳細(xì)分析了有限差分的時(shí)間和空間頻散特性及其對(duì)成像和反演精度的影響；發(fā)展了抗時(shí)間頻散、抗空間頻散、抗時(shí)間+空間頻散的逆時(shí)偏移和全波形反演方法；通過理論數(shù)據(jù)和實(shí)際資料對(duì)提出的成像和反演方法進(jìn)行了測(cè)試.

1 有限差分的時(shí)間和空間頻散

二階常密度聲波方程形式如下(Dablain，1986；Liu and Sen，2009b)：

(1)

式中，L=?2/?x2+?2/?z2，v為傳播速度，p為壓強(qiáng)，s為震源項(xiàng).常規(guī)有限差分時(shí)間和空間偏導(dǎo)數(shù)離散格式為：

(2)

(3)

其中，cn為差分系數(shù)，M為空間算子長(zhǎng)度參數(shù)，Δt為時(shí)間步長(zhǎng)，Δx和Δz為空間采樣間隔.方程(2)和(3)兩邊進(jìn)行傅里葉變換并化解可得：

(4)

(5)

其中，ω為角頻率，k為波數(shù)，kx和kz分別為水平與垂直方向的波數(shù).

(6)

(7)

圖1 常規(guī)有限差分法離散前后的頻率(a)和波數(shù)(b)變化曲線Fig.1 Variation curves of frequency (a) and wavenumber (b) before and after the discretization in the conventional finite difference method

(8)

有限差分法的空間精度隨M的增加而增大.當(dāng)M趨向于無窮大時(shí)，有限差分法變?yōu)閭巫V法，空間偏導(dǎo)數(shù)通過式(9)求解(Kosloff and Baysal，1982)：

Lp=F-1(F(Lp)),

(9)

其中，F(xiàn)和F-1分別表示二維正反傅里葉變換.與空間頻散不同，時(shí)間頻散是可預(yù)測(cè)的，與傳播路徑和速度等因素?zé)o關(guān)(Koene et al.，2018).有限差分的時(shí)間頻散是由于數(shù)值離散引起的頻率移動(dòng)(從ω到ωFD)而產(chǎn)生的.將具有頻散的數(shù)據(jù)從ωFD?；氐秸鎸?shí)頻率ω處即可實(shí)現(xiàn)時(shí)間頻散校正，得到原始的無頻散數(shù)據(jù).由方程(6)可得：

(10)

方程(10)定量表示離散前后的頻率變化關(guān)系.時(shí)間頻散校正步驟為(Wang and Xu，2015)：

(1)通過方程(10)計(jì)算不同ω對(duì)應(yīng)的ωFD；

常規(guī)有限差分具有時(shí)間2階和空間2M階精度，同時(shí)產(chǎn)生時(shí)間和空間頻散.偽譜法的空間模擬精度為無窮階.采用時(shí)間有限差分和偽譜法(方程(2)和(9))求解波動(dòng)方程，可以得到只包含時(shí)間頻散的數(shù)據(jù).時(shí)間頻散校正法能夠消除常規(guī)時(shí)間二階有限差分引入的時(shí)間頻散.采用時(shí)間頻散校正和空間有限差分(方程(3)和(10))可以得到只包含空間頻散的數(shù)據(jù).將時(shí)間頻散校正和偽譜法(方程(9)和(10))相結(jié)合可獲取不含時(shí)間和空間頻散的數(shù)據(jù).為了書寫方便，下文中將有限差分、偽譜法和時(shí)間頻散校正分別記為FD、PS和TDC.因此，F(xiàn)D、PS、FD+TDC和PS+TDC分別表示全頻散、只有時(shí)間頻散、只有空間頻散和無頻散數(shù)據(jù)或者對(duì)應(yīng)的方法.

這里采用一個(gè)均勻介質(zhì)模型研究有限差分的時(shí)間和空間頻散特性.介質(zhì)速度為2000 m·s-1，網(wǎng)格點(diǎn)數(shù)為301×301，空間網(wǎng)格大小為10 m×10 m，時(shí)間步長(zhǎng)為 2.0 ms，最大記錄時(shí)間為2.0 s.震源為30 Hz的雷克子波，位于點(diǎn)(1000 m,1000 m)處.圖2給出不同方法得到的點(diǎn)(3000 m,2000 m)處的波形圖及頻譜.其中，Semi-analytical為由Cagniard-De Hoop方法(De Hoop，1960)得到的半解析波場(chǎng).FD表示時(shí)間2階空間10階精度的有限差分.由圖2可見：時(shí)間頻散(PS)會(huì)引起波場(chǎng)相移，空間頻散(FD+TDC)導(dǎo)致波形震蕩；時(shí)間和空間頻散分別出現(xiàn)在正常波前的前部和尾部；PS+TDC可以有效壓制時(shí)間和空間頻散，獲得與解析解相近的模擬結(jié)果；空間頻散引起頻譜畸變，時(shí)間頻散對(duì)頻譜的影響相對(duì)較小.

圖2 均勻介質(zhì)模型(3000 m, 2000 m)處不同方法的波形圖(a)及頻譜(b)Fig.2 Waveforms (a) and spectra (b) at the point (3000 m, 2000 m) by different methods on a homogeneous model

2 抗頻散逆時(shí)偏移方法

逆時(shí)偏移包括震源波場(chǎng)正向傳播、檢波點(diǎn)波場(chǎng)反向傳播和正反向波場(chǎng)互相關(guān)成像.將FD、PS、FD+TDC和PS+TDC應(yīng)用于波場(chǎng)正傳和反傳中，即可研究時(shí)間和空間頻散對(duì)逆時(shí)偏移的影響.但直接進(jìn)行時(shí)間頻散校正，需要存儲(chǔ)所有網(wǎng)格點(diǎn)所有時(shí)刻的正反向波場(chǎng)，并進(jìn)行大量的正反傅里葉變換，計(jì)算和存儲(chǔ)需求巨大.為了提高計(jì)算效率，Koene等(2018)介紹了一種實(shí)用的抗時(shí)間頻散逆時(shí)偏移(FD+TDC)方法：震源s中加入時(shí)間頻散，源波場(chǎng)正向傳播；觀測(cè)波場(chǎng)中加入時(shí)間頻散，檢波點(diǎn)波場(chǎng)反向傳播；正反向波場(chǎng)互相關(guān)成像.向無頻散數(shù)據(jù)中添加時(shí)間頻散就是實(shí)現(xiàn)從ω到ωFD的過程，具體步驟為(Koene et al.，2018)：

(1)通過方程(6)計(jì)算不同ωFD對(duì)應(yīng)的ω；

震源波場(chǎng)中只包含源-反射點(diǎn)路徑上的時(shí)間頻散.觀測(cè)波場(chǎng)中加入時(shí)間頻散后包含源-反射點(diǎn)-檢波點(diǎn)整個(gè)路徑的時(shí)間頻散.由于波場(chǎng)正向和反向傳播沿不同的時(shí)間方向積累時(shí)間頻散，檢波點(diǎn)波場(chǎng)反向傳播會(huì)抵消掉反射點(diǎn)-檢波點(diǎn)路徑上積累的時(shí)間頻散.因此，震源波場(chǎng)和檢波點(diǎn)波場(chǎng)具有相同的時(shí)間頻散(相移)，仍可以實(shí)現(xiàn)零延遲互相關(guān)成像.

偽譜法具有空間無窮階精度，采用偽譜法進(jìn)行波場(chǎng)正反向延拓得到抗空間頻散逆時(shí)偏移(PS)方法.將時(shí)間頻散校正和偽譜法結(jié)合得到抗時(shí)間+空間頻散逆時(shí)偏移(PS+TDC)方法：震源s中加入時(shí)間頻散，偽譜法源波場(chǎng)正向傳播；觀測(cè)波場(chǎng)中加入時(shí)間頻散，偽譜法檢波點(diǎn)波場(chǎng)反向傳播；正反向波場(chǎng)互相關(guān)成像.

3 抗頻散全波形反演方法

全波形反演包括震源波場(chǎng)正向傳播、殘差波場(chǎng)反向傳播、模型參數(shù)梯度計(jì)算和迭代更新.將FD、PS、FD+TDC和PS+TDC應(yīng)用于波場(chǎng)正傳和反傳中，同樣可以研究時(shí)間和空間頻散對(duì)全波形反演的影響.抗時(shí)間頻散全波形反演(FD+TDC)方法為：震源s中加入時(shí)間頻散，源波場(chǎng)正向傳播；觀測(cè)波場(chǎng)中加入時(shí)間頻散，殘差波場(chǎng)反向傳播；模型參數(shù)梯度計(jì)算；反演迭代(如共軛梯度法)直至滿足收斂條件.

由于震源波場(chǎng)和伴隨波場(chǎng)中殘留有源-反射點(diǎn)路徑的時(shí)間頻散，導(dǎo)致得到的模型參數(shù)梯度與無頻散情況仍有差異.模擬波場(chǎng)和觀測(cè)波場(chǎng)包含源-反射點(diǎn)-檢波點(diǎn)路徑的時(shí)間頻散引起波形震蕩、頻譜改變，也會(huì)增加反演的多解性.因此，抗空間頻散(PS)和抗時(shí)間+空間頻散(PS+TDC)全波形反演方法仍無法徹底消除時(shí)間頻散的影響.與波場(chǎng)延拓相比，觀測(cè)波場(chǎng)中添加時(shí)間頻散的計(jì)算量可忽略不計(jì)(且可在偏移和反演之前完成).抗時(shí)間+空間頻散(PS+TDC)逆時(shí)偏移和全波形反演方法的計(jì)算效率與空間頻散校正(PS)方法相近.

4 數(shù)值算例

采用salt模型、Sigsbee2A模型、實(shí)際資料對(duì)抗頻散逆時(shí)偏移和全波形反演方法進(jìn)行測(cè)試.通過對(duì)比抗空間頻散和抗時(shí)間+空間頻散方法的結(jié)果研究時(shí)間頻散對(duì)成像和反演的影響；對(duì)比抗時(shí)間頻散和抗時(shí)間+空間頻散方法的結(jié)果研究空間頻散對(duì)成像和反演的影響.基于不同方法與參考解的接近程度驗(yàn)證抗時(shí)間+空間頻散方法的有效性.

4.1 逆時(shí)偏移

SEG/EAGEsalt模型如圖3所示.空間網(wǎng)格大小為20 m×20 m，時(shí)間步長(zhǎng)為2.0 ms，最大記錄時(shí)間為4.0 s.震源采用主頻為30 Hz的Ricker子波.31個(gè)震源和601個(gè)檢波點(diǎn)均勻分布于地表.采用偽譜法和小時(shí)間步長(zhǎng)(Δt=0.5 ms)模擬無頻散觀測(cè)記錄.FD表示時(shí)間2階空間16階精度的有限差分.常規(guī)(FD)和抗空間頻散(PS)逆時(shí)偏移中直接對(duì)無頻散觀測(cè)記錄進(jìn)行反傳，而抗時(shí)間頻散(FD+TDC)和抗時(shí)間+空間頻散(PS+TDC)方法中先在觀測(cè)記錄中加入時(shí)間頻散后再反傳.圖4為第16炮(7800 m, 2160 m)處不同偏移方法的正反向傳播波場(chǎng)和頻譜.為便于比較，采用偽譜法和小時(shí)間步長(zhǎng)(Δt=0.5 ms)進(jìn)行波場(chǎng)正反向延拓并成像作為參考解(Ref).由圖4可見，正傳和反傳波場(chǎng)的時(shí)間頻散分別引起波形向負(fù)時(shí)間和正時(shí)間方向移動(dòng)，正傳和反傳波場(chǎng)的空間頻散分別出現(xiàn)在正常波前時(shí)間增大和減小的方向；常規(guī)(FD)和抗時(shí)間頻散(FD+TDC)方法的正傳波場(chǎng)同時(shí)包含時(shí)間和空間頻散(相移和震蕩)，而抗空間頻散(PS)和抗時(shí)間+空間頻散(PS+TDC)方法的正傳波場(chǎng)中只積累了源-反射點(diǎn)路徑的時(shí)間頻散(相移)；常規(guī)(FD)和抗空間頻散(PS)方法的反傳波場(chǎng)包含檢波點(diǎn)-反射點(diǎn)路徑的時(shí)間頻散，而抗時(shí)間頻散(FD+TDC)和抗時(shí)間+空間頻散(PS+TDC)方法的反傳波場(chǎng)中積累了源-反射點(diǎn)路徑的時(shí)間頻散，常規(guī)(FD)和抗時(shí)間頻散(FD+TDC)方法的反傳波場(chǎng)空間頻散較強(qiáng)；常規(guī)(FD)和抗時(shí)間頻散(FD+TDC)方法的正傳波場(chǎng)頻譜變形嚴(yán)重，抗空間頻散(PS)和抗時(shí)間+空間頻散(PS+TDC)方法的正反傳波場(chǎng)頻譜與參考解的差異較小.圖5給出不同方法的正反向傳播波場(chǎng)匹配圖.常規(guī)(FD)和抗空間頻散(PS)方法正反傳波場(chǎng)中分別積累了源-反射點(diǎn)和檢波點(diǎn)-反射點(diǎn)路徑的時(shí)間頻散，波形沿相反的方向移動(dòng)，導(dǎo)致正反傳波場(chǎng)無法匹配；抗時(shí)間頻散(FD+TDC)和抗時(shí)間+空間頻散(PS+TDC)方法的正反傳波場(chǎng)中包含相同的源-反射點(diǎn)路徑時(shí)間頻散，波形都沿時(shí)間負(fù)方向移動(dòng)，正反傳波場(chǎng)相位相同，仍然能進(jìn)行互相關(guān)條件成像；采用偽譜法(PS)進(jìn)行正反向波場(chǎng)延拓，可以有效壓制空間頻散，減小波形震蕩對(duì)成像精度的影響.

圖3 鹽丘模型Fig.3 Salt dome model

圖4 鹽丘模型(7800 m, 2160 m)處不同方法的正反向波場(chǎng)及頻譜(a) 正向波場(chǎng)； (b) 反向波場(chǎng)； (c) 正向波場(chǎng)頻譜； (d) 反向波場(chǎng)頻譜.Fig.4 Wavefields and spectra at the point (7800 m, 2160 m) by different methods on salt dome model(a) Forward wavefield； (b) Backward wavefield； (c) Spectrum of (a)； (d) Spectrum of (b).

圖5 鹽丘模型不同方法的正反向波場(chǎng)匹配圖(a) Ref; (b) FD; (c) PS; (d) FD+TDC; (e) PS+TDC.Fig.5 Matching diagrams of forward- and backward-propagated wavefields by different methods for the salt dome model

圖6為不同偏移方法得到的成像剖面.由圖6可知，常規(guī)(FD)方法同時(shí)受時(shí)間和空間頻散影響，無法準(zhǔn)確成像；抗空間頻散(PS)方法只與時(shí)間頻散有關(guān)，出現(xiàn)同相軸不聚焦、分辨率降低、反射界面偏離真實(shí)位置等現(xiàn)象；抗時(shí)間頻散(FD+TDC)方法只受空間頻散影響，產(chǎn)生了許多高頻噪聲(如黑色箭頭所示)和虛假反射界面(如白色箭頭所示)；抗時(shí)間+空間頻散(PS+TDC)方法幾乎不受頻散影響，得到了與參考解相近的結(jié)果.圖7給出不同方法像的波數(shù)譜.抗空間頻散(PS)方法的譜中包含垂直方向的線條(如黑色箭頭所示)，對(duì)應(yīng)于空間域的周期性虛假水平同相軸；抗時(shí)間頻散(FD+TDC)方法的譜中出現(xiàn)了一些高波數(shù)成分(如黑色箭頭所示)，對(duì)應(yīng)于空間域的高頻成像噪聲；抗時(shí)間頻散(FD+TDC)方法的波數(shù)譜與參考解能夠較好地吻合.

圖6 鹽丘模型不同方法的成像結(jié)果(a) Ref; (b) FD; (c) PS; (d) FD+TDC; (e) PS+TDC.Fig.6 Images by different methods for the salt dome model

圖7 鹽丘模型不同方法成像結(jié)果的波數(shù)譜(a) Ref; (b) FD; (c) PS; (d) FD+TDC; (e) PS+TDC.Fig.7 Spectra of images by different methods for the salt dome model

下面采用海上某工區(qū)實(shí)際地震資料對(duì)不同的偏移方法進(jìn)行測(cè)試.該數(shù)據(jù)包含240炮，每炮120道，炮間距50 m，道間距25 m.網(wǎng)格大小為25 m×25 m，時(shí)間步長(zhǎng)為2.0 ms，最大記錄時(shí)間為4.0 s.圖8給出其中5炮的地震記錄.圖9為偏移速度模型.震源子波是通過軟件提取的最小相位子波.FD表示時(shí)間2階空間10階精度的有限差分.圖10展示了不同偏移方法的成像結(jié)果.由圖10可見：常規(guī)(FD)和抗時(shí)間頻散(FD+TDC)方法偏移剖面中存在較強(qiáng)的成像噪聲；抗空間頻散(PS)方法偏移剖面同相軸的連續(xù)性較差；抗時(shí)間+空間頻散(FD+TDC)方法不受時(shí)間和空間頻散的影響，得到同相軸連續(xù)性更好、噪聲少、分辨率更高的像.圖11給出不同方法成像結(jié)果的波數(shù)譜.通過時(shí)間頻散校正和偽譜法壓制時(shí)間和空間頻散，可有效恢復(fù)出缺失的波數(shù)成分并削弱成像噪聲(如虛線框和箭頭所示).但實(shí)際應(yīng)用中速度模型和震源子波未知，成像結(jié)果同時(shí)受速度和震源精度以及數(shù)值頻散的影響，導(dǎo)致不同方法的偏移效果差異不明顯.只有消除其他不利因素(如速度誤差、震源誤差、照明不均等)影響，本文發(fā)展的抗頻散逆時(shí)偏移方法才更加有效.

圖8 實(shí)際地震記錄Fig.8 Real seismic record

圖9 實(shí)際資料偏移速度模型Fig.9 Migration velocity model for real data

圖10 實(shí)際資料不同方法的成像結(jié)果(a) FD; (b) PS; (c FD+TDC; (d) PS+TDC.Fig.10 Images by different methods for real data

圖11 實(shí)際資料不同方法的成像結(jié)果波數(shù)譜(a) FD; (b) PS; (c) FD+TDC; (d) PS+TDC.Fig.11 Spectra of images by different methods for real data

4.2 全波形反演

Sigsbee2A模型如圖12所示.計(jì)算區(qū)域?yàn)?600 m×2960 m，網(wǎng)格大小為16 m×16 m，時(shí)間步長(zhǎng)為1.5 ms，最大記錄長(zhǎng)度為3.0 s.震源采用主頻為30 Hz的Ricker子波.36個(gè)震源和351個(gè)檢波點(diǎn)均勻分布于地表.初始速度模型為真實(shí)模型的平滑結(jié)果，如圖13所示.由于初始模型與真實(shí)模型相近，我們進(jìn)行單尺度反演，不同方法迭代10次.采用偽譜法和小時(shí)間步長(zhǎng)(Δt=0.375 ms)模擬無頻散觀測(cè)記錄.抗時(shí)間頻散(FD+TDC)和抗時(shí)間+空間頻散(PS+TDC)方法中先在觀測(cè)記錄中加入時(shí)間頻散再計(jì)算殘差并進(jìn)行波場(chǎng)反向延拓.圖14為第18炮(4000 m,2240 m)處不同方法第一次迭代時(shí)的正反向傳播波場(chǎng)和頻譜.FD表示時(shí)間2階空間12階精度的有限差分.采用偽譜法和小時(shí)間步長(zhǎng)(Δt=0.375 ms)進(jìn)行波場(chǎng)正反向延拓并反演作為參考解(Ref).由圖14可見：常規(guī)(FD)、抗空間頻散(PS)、抗時(shí)間頻散(FD+TDC)和抗時(shí)間+空間頻散(PS+TDC)方法的正傳波場(chǎng)相對(duì)于參考解向時(shí)間負(fù)方向發(fā)生相移；常規(guī)(FD)和抗空間頻散(PS)的反傳波場(chǎng)相對(duì)于參考解向時(shí)間正方向發(fā)生相移，而抗時(shí)間頻散(FD+TDC)和抗時(shí)間+空間頻散PS+TDC)方法的反傳波場(chǎng)仍向時(shí)間負(fù)方向移動(dòng)；正反向波場(chǎng)的空間頻散分別出現(xiàn)在正常波前時(shí)間增大和減小的方向，偽譜法(PS)可以有效壓制空間頻散.圖15給出不同反演方法第一次迭代時(shí)的正反傳波場(chǎng)匹配圖.時(shí)間頻散校正可以保證正反向波場(chǎng)中包含等量的(源-反射點(diǎn)路徑)時(shí)間頻散.抗時(shí)間頻散(FD+TDC)和抗時(shí)間+空間頻散(PS+TDC)方法的正反傳波場(chǎng)的時(shí)移量與參考解正反傳波場(chǎng)的時(shí)移量大致相同.正反傳波場(chǎng)的不匹配主要是由于速度存在誤差導(dǎo)致的.

圖12 Sigsbee2A模型Fig.12 Sigsbee2A model

圖13 Sigsbee2A模型初始模型Fig.13 Initial model for the Sigsbee2A model

圖14 Sigsbee2A模型(4000 m, 2240 m)處不同方法的正反向波場(chǎng)及頻譜(a) 正向波場(chǎng)； (b) 反向波場(chǎng)； (c) 正向波場(chǎng)頻譜； (d) 反向波場(chǎng)頻譜.Fig.14 Forward and backward wavefields and frequency spectra at the point (4000 m, 2240 m) by different methods for the Sigsbee2A model(a) Forward wavefield； (b) Backward wavefield； (c) Spectrum of (a)； (d)Spectrum of (b).

圖15 Sigsbee2A模型不同方法的正反向波場(chǎng)匹配圖(a) Ref; (b) FD; (c) PS; (d) FD+TDC; (e) PS+TDC.Fig.15 Matching diagrams of forward- and backward-propagated wavefields by different methods for the Sigsbee2A model

圖16為不同反演方法得到的速度模型.由圖16可知，常規(guī)(FD)和抗時(shí)間頻散(FD+TDC)方法受空間頻散影響嚴(yán)重，無法得到滿意的結(jié)果；抗空間頻散(PS)和抗時(shí)間+空間頻散(PS+TDC)方法反演的速度模型有一定改善.為進(jìn)一步比較，圖17給出Distance = 2720 m處不同方法的反演結(jié)果.抗空間頻散(PS)方法受時(shí)間頻散影響，反演精度較低；抗時(shí)間+空間頻散(PS+TDC)方法可以緩解時(shí)間和空間頻散對(duì)反演的影響，獲得了與參考解更接近的結(jié)果.圖18展示了不同方法的收斂曲線.強(qiáng)的空間頻散導(dǎo)致常規(guī)(FD)和抗時(shí)間頻散(FD+TDC)方法無法收斂；抗時(shí)間+空間頻散(PS+TDC)方法受時(shí)間頻散影響較小，比抗空間頻散(PS)方法的收斂速度快.抗時(shí)間+空間頻散(PS+TDC)方法中震源和伴隨波場(chǎng)中殘留有源-反射點(diǎn)路徑的時(shí)間頻散，模擬和觀測(cè)波場(chǎng)也包含源-反射點(diǎn)-檢波點(diǎn)路徑的時(shí)間頻散.因此，抗頻散(PS+TDC)全波形反演方法仍無法徹底消除時(shí)間頻散的影響，其反演結(jié)果和收斂性與參考解仍有一定距離.

圖16 Sigsbee2A模型不同方法的反演結(jié)果(a) Ref; (b) FD; (c) PS; (d) FD+TDC; (e) PS+TDC.Fig.16 Inversion results by different methods for the Sigsbee2A model

圖17 Sigsbee2A模型Distance=2720 m處不同方法的反演結(jié)果Fig.17 Inversion results at distance 2720 m by different methods for the Sigsbee2A model

圖18 Sigsbee2A模型不同方法的收斂曲線Fig.18 Convergence curves by different methods for the Sigsbee2A model

下面采用一維初始模型(如圖19所示)對(duì)不同的反演方法進(jìn)行測(cè)試.由于初始模擬與真實(shí)模型相差較遠(yuǎn)，我們進(jìn)行多尺度反演：0～5 Hz、0～10 Hz和0～15 Hz，每個(gè)尺度迭代10次.圖20給出Distance=2720 m處不同方法的反演結(jié)果.在大尺度低頻反演中，時(shí)間和空間頻散相對(duì)較弱，對(duì)全波形反演的影響非常??；通過時(shí)間頻散校正(TDC)減小時(shí)間頻散的影響，可以微弱改善反演精度.

圖19 Sigsbee2A模型一維初始模型Fig.19 1D initial model for the Sigsbee2A model

圖20 一維初始模型時(shí)Sigsbee2A模型Distance=2720 m處不同方法的反演結(jié)果Fig.20 Inversion results at distance 2720 m by different methods with the 1D initial model for the Sigsbee2A model

5 結(jié)論

本文詳細(xì)研究了有限差分的時(shí)間和空間頻散特性及其對(duì)逆時(shí)偏移和全波形反演的影響.通過時(shí)間有限差分+偽譜法和時(shí)間頻散校正+空間有限差分分離出時(shí)間和空間頻散，進(jìn)而發(fā)展了抗時(shí)間頻散、抗空間頻散、抗時(shí)間+空間頻散的逆時(shí)偏移和全波形反演方法.模型算例和實(shí)際資料應(yīng)用表明：逆時(shí)偏移同時(shí)受時(shí)間和空間頻散影響，時(shí)間頻散導(dǎo)致同相軸連續(xù)性變差、成像位置偏離，空間頻散產(chǎn)生高頻噪聲和虛假反射界面；全波形反演在低頻大尺度反演中幾乎不受時(shí)間和空間頻散影響，高頻精細(xì)反演中時(shí)間頻散會(huì)降低反演精度，空間頻散可能導(dǎo)致不收斂；抗頻散方法可以有效緩解時(shí)間和空間頻散影響，獲得滿意的偏移剖面和反演結(jié)果.

時(shí)間頻散校正和偽譜法不依賴于差分網(wǎng)格和波動(dòng)方程，本文的抗頻散方法可以直接推廣到復(fù)雜介質(zhì)(如變密度聲波、彈性、各向異性等)逆時(shí)偏移和全波形反演中.當(dāng)模型復(fù)雜且存在速度劇烈變化界面時(shí)，偽譜法的模擬精度變差.通過對(duì)速度進(jìn)行平滑或者采用譜元法可以緩解上述問題.時(shí)間頻散校正+優(yōu)化空間有限差分也是一種兼顧精度和計(jì)算效率的方法.為徹底消除時(shí)間頻散對(duì)全波形反演的影響，需要進(jìn)一步研究基于高階時(shí)間有限差分的抗頻散方法.

Dong L G, Ma Z T, Cao J Z, et al. 2000. A staggered-grid high-order difference method of one-order elastic wave equation.ChineseJournalofGeophysics(in Chinese), 43(3): 411-419.

Liu X, Liu Y, Ren Z M, et al. 2019. Perfectly matched layer boundary conditions for frequency-domain acoustic wave simulation in the mesh-free discretization.GeophysicalProspecting, 67(7): 1732-1744.

Tan S R, Huang L J. 2014. An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation.GeophysicalJournalInternational, 197(2): 1250-1267.