近似常Q模型的黏聲各向異性純qP波逆時(shí)偏移

張閃閃， 谷丙洛， 丁奕文， 王霽川， 李振春

中國(guó)石油大學(xué)(華東)深層油氣重點(diǎn)實(shí)驗(yàn)室， 青島 266580

0 引言

地下真實(shí)介質(zhì)中普遍存在著各向異性(Thomsen，1986)，地震波在各向異性介質(zhì)中傳播時(shí)，其彈性特征會(huì)隨著傳播方向發(fā)生變化，具體表現(xiàn)為地震波速度隨傳播方向變化而變化(Xu et al.，2020).因此利用各向同性逆時(shí)偏移方法對(duì)各向異性地震資料進(jìn)行成像時(shí)，由于地震波走時(shí)不準(zhǔn)確，會(huì)導(dǎo)致構(gòu)造位置成像不準(zhǔn)、同相軸能量不聚焦、繞射波收斂不完全等問題的出現(xiàn).此外，理論和實(shí)驗(yàn)證明實(shí)際地下介質(zhì)是黏彈性的(Robertsson et al.，1994).地震波在黏彈性介質(zhì)中傳播時(shí)，地下介質(zhì)固有的黏滯效應(yīng)會(huì)造成地震波能量衰減和相位畸變.在強(qiáng)衰減介質(zhì)中，受介質(zhì)黏滯性影響，實(shí)際采集的地震資料將存在明顯的振幅衰減和相位畸變，采用無吸收衰減補(bǔ)償?shù)钠瞥上穹椒ㄟM(jìn)行成像時(shí)，偏移剖面分辨率將大幅降低——尤其在強(qiáng)吸收介質(zhì)下方的地層，導(dǎo)致油氣預(yù)測(cè)出現(xiàn)較大誤差(Aki and Richards，2002；Shen et al.，2018).

Alkhalifah(1998)提出各向異性聲學(xué)近似假設(shè)并推導(dǎo)了適用于VTI介質(zhì)的四階偽聲波方程，但該方程包含波場(chǎng)的四階偏導(dǎo)數(shù)，數(shù)值求解對(duì)內(nèi)存要求很高.為了降低計(jì)算成本，Du等(2008)利用輔助算子將偽聲波方程簡(jiǎn)化為兩個(gè)耦合的二階偏微分方程，實(shí)現(xiàn)qP波的高效數(shù)值模擬.然而，在求解基于聲學(xué)近似的TTI介質(zhì)波動(dòng)方程時(shí)，由于將對(duì)稱軸上的qSV波速度設(shè)置為零，導(dǎo)致求解過程存在數(shù)值不穩(wěn)定(Grechka et al.，2004).針對(duì)此問題，F(xiàn)letcher等(2009)假設(shè)橫波速度大于零，引入與全彈性頻散關(guān)系一致的二階項(xiàng)來修正Alkhalifah(1998)所提出的頻散關(guān)系，并通過坐標(biāo)變換將VTI介質(zhì)qP波方程推廣到TTI介質(zhì).Zhang等(2011)在聲學(xué)近似條件下，使用伴隨算子和協(xié)變導(dǎo)數(shù)算子將波場(chǎng)變量旋轉(zhuǎn)到局部?jī)A斜坐標(biāo)系中，直接利用VTI介質(zhì)中的彈性波方程推導(dǎo)出穩(wěn)定的TTI波動(dòng)方程.為了消除震源處的偽橫波干擾，Yoon等(2010)提出了橢圓各向異性近似方法，即在震源周圍設(shè)置“ε=δ”的光滑圓錐形區(qū)域，該方法可有效抑制偽橫波假象.然而，這種方法雖然旨在減少橫波偽影，卻改變了地震波的運(yùn)動(dòng)學(xué)特征，并在縱波模擬中保留了人工橫波分量.為了避免這些現(xiàn)象，Zhang(2009)通過濾波投影方法來壓制偽橫波，以提高成像剖面精度.程玖兵等(2013)推導(dǎo)出適用于任意各向異性介質(zhì)的qP波偽純模式波動(dòng)方程，并利用偏振矢量將縱、橫波分離，進(jìn)一步降低計(jì)算成本、壓制橫波干擾.Xu和Zhou(2014)基于聲學(xué)近似假設(shè)，提出了VTI介質(zhì)qP波方程，該方程利用偽差分算子將qP波和qS波解耦表示，其數(shù)值模擬結(jié)果不含偽橫波假象.劉文卿等(2016)推導(dǎo)出可用顯式有限差分格式求解的TTI介質(zhì)下qP波控制方程，并提出正則化有限橫波方法解決波場(chǎng)傳播及數(shù)值計(jì)算的不穩(wěn)定問題.黃建平等(2016)推導(dǎo)了TTI介質(zhì)一階速度-應(yīng)力方程并用高階有限差分格式求解，同時(shí)引入了基于PML原理的一階qP波方程的偽譜法遞推格式來實(shí)現(xiàn)波場(chǎng)的穩(wěn)定延拓.Cheng和Kang(2016)等從物理學(xué)角度出發(fā)，利用波矢量方向與偽橫波震動(dòng)方向的偏差，從耦合的波場(chǎng)中提取出標(biāo)量純qP波場(chǎng)，以提高逆時(shí)偏移成像精度.郭成峰等(2017)通過對(duì)Harlan公式中的非橢圓項(xiàng)進(jìn)行修正推導(dǎo)出TTI介質(zhì)純P波方程，并利用偽譜法和有限差分法聯(lián)合實(shí)現(xiàn)高精度波場(chǎng)延拓及逆時(shí)偏移.楊鵬等(2017)推導(dǎo)出解耦的TTI介質(zhì)純qP波近似方程，將其中的偏微分算子分解成一個(gè)拉普拉斯算子和一個(gè)標(biāo)量算子，其數(shù)值模擬結(jié)果可有效實(shí)現(xiàn)振幅均衡.郭旭等(2019)利用Poynting矢量分離地震方向行波，并在復(fù)雜VTI介質(zhì)中成像，該方法在偏移速度不準(zhǔn)確時(shí)仍能有效成像.何兵壽等(2021)指出，在TI介質(zhì)中的qP波逆時(shí)偏移領(lǐng)域需要對(duì)qP波方程的動(dòng)力學(xué)精度，qP波逆時(shí)偏移中的隨機(jī)邊界技術(shù)及波場(chǎng)傳播方向的求取方法深入研究.

除了各向異性之外，學(xué)者們?cè)谒p補(bǔ)償方面也做出了許多努力.Kjartansson(1979)提出一種不依賴頻率的Q值模型，又稱為常Q模型.基于常Q模型，Carcione等(2002)推導(dǎo)出包含分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)的黏聲波動(dòng)方程來描述聲波在衰減介質(zhì)下的傳播，但該方程需要存儲(chǔ)所有歷史時(shí)刻的波場(chǎng)信息，要求極大的存儲(chǔ)空間.在此基礎(chǔ)上，Chen和Holm(2004)將分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)轉(zhuǎn)化為分?jǐn)?shù)階拉普拉斯算子，通過偽譜法在波數(shù)域進(jìn)行數(shù)值求解，提高了計(jì)算效率，節(jié)約了內(nèi)存空間.Treeby和Cox(2010)將振幅衰減項(xiàng)和相位頻散項(xiàng)解耦表示，方便在成像過程中進(jìn)行衰減補(bǔ)償.基于近似常Q模型，Yan和Liu(2013)采用優(yōu)化的時(shí)空域高階有限差分算法對(duì)黏聲波動(dòng)方程進(jìn)行數(shù)值求解，實(shí)現(xiàn)了高精度疊前逆時(shí)偏移.Zhu等(2014)提出含有分?jǐn)?shù)階拉普拉斯算子的時(shí)域常Q黏聲波方程，將振幅衰減項(xiàng)和相位頻散項(xiàng)解耦表示，進(jìn)行穩(wěn)定補(bǔ)償?shù)酿ぢ暷鏁r(shí)偏移.Li 等(2016)提出利用頻散關(guān)系將常Q黏聲波方程中的振幅衰減和相位畸變進(jìn)一步分離，使其能夠應(yīng)用于Q值較小地層，有效改善成像效果.吳玉等(2017)基于解耦的分?jǐn)?shù)階拉普拉斯算子黏聲波方程，采用交錯(cuò)網(wǎng)格有限差分法實(shí)現(xiàn)數(shù)值模擬和逆時(shí)偏移，在波場(chǎng)逆時(shí)延拓過程中同時(shí)完成頻散校正和衰減補(bǔ)償，大幅提高計(jì)算效率和成像質(zhì)量.周彤等(2018)采用多項(xiàng)式多級(jí)優(yōu)化方法校正衰減引起的相位畸變，提高逆時(shí)偏移成像精度.在時(shí)間域二階位移形式的常分?jǐn)?shù)階拉普拉斯算子黏聲波方程的基礎(chǔ)上，陳漢明等(2020)推導(dǎo)了一階速度-應(yīng)力形式的常分?jǐn)?shù)階拉普拉斯算子黏聲波動(dòng)方程，并提出一種適用于該方程的卷積型完全匹配層吸收邊界加載方法.

在黏聲各向同性波動(dòng)方程的基礎(chǔ)上，學(xué)者們逐漸將其向各向異性介質(zhì)發(fā)展.Suh等(2012)將Zhang等(2010)推導(dǎo)的含有偽差分算子的黏聲波動(dòng)方程擴(kuò)展到VTI介質(zhì)，但該方程忽略了衰減引起的相位頻散，只對(duì)振幅衰減進(jìn)行補(bǔ)償.嚴(yán)紅勇和劉洋(2012)以Carcione黏彈各向異性理論為基礎(chǔ)，推導(dǎo)了適用于黏彈TTI介質(zhì)的二維三分量一階速度-應(yīng)力方程，同時(shí)采用旋轉(zhuǎn)交錯(cuò)網(wǎng)格任意偶數(shù)階精度有限差分格式進(jìn)行求解，實(shí)現(xiàn)了該類介質(zhì)的地震波場(chǎng)數(shù)值模擬.Xie等(2015)根據(jù)線性黏彈介質(zhì)的應(yīng)力-應(yīng)變關(guān)系推導(dǎo)出TTI介質(zhì)中的黏聲波動(dòng)方程，并發(fā)展了TTI介質(zhì)中的Q補(bǔ)償成像方法，有效改善剖面質(zhì)量.Zhu(2017)推導(dǎo)了適用于任意各向異性速度和衰減的時(shí)域位移-應(yīng)力黏彈波動(dòng)方程，該方程使用分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)表征衰減各向異性，因此需要存儲(chǔ)所有歷史時(shí)刻的波場(chǎng)信息，存儲(chǔ)量較大.為了降低計(jì)算成本，Hu等(2017)發(fā)展了基于有限差分解法的TI介質(zhì)中的Q補(bǔ)償逆時(shí)偏移成像方法，該方法大幅提高了計(jì)算效率，便于代碼并行化和GPU實(shí)現(xiàn).Bai等(2019)提出VTI介質(zhì)中振幅衰減和相位頻散解耦的Q補(bǔ)償逆時(shí)偏移成像算法，通過改變振幅衰減項(xiàng)前的符號(hào)、保持相位頻散項(xiàng)不變來進(jìn)行衰減補(bǔ)償，大幅提高VTI介質(zhì)的成像質(zhì)量.在弱衰減假設(shè)條件下，Zhu和Bai(2019)提出了振幅衰減和相位頻散解耦的分?jǐn)?shù)階拉普拉斯黏聲VTI波動(dòng)方程，該方程在保持模擬精度的同時(shí)可以大幅提高數(shù)值模擬的效率.Hao和Alkhalifah(2019)對(duì)正交各向異性介質(zhì)中標(biāo)量和矢量黏聲波動(dòng)方程的一般表示形式進(jìn)行了推導(dǎo)，使其在多個(gè)衰減模型下均可靈活適用.Zhou等(2020)將線性黏滯本構(gòu)關(guān)系擴(kuò)展為一系列整數(shù)階微分項(xiàng)和一個(gè)唯一的積分項(xiàng)，實(shí)現(xiàn)了振幅和相位的解耦，并提出基于有限差分解法的Q補(bǔ)償逆時(shí)偏移算法框架，該成像方法可在更寬的頻率范圍內(nèi)實(shí)現(xiàn)精確的衰減補(bǔ)償.程時(shí)俊等(2021)將衰減各向異性和速度各向異性與常Q模型相結(jié)合，建立了黏彈性衰減VTI介質(zhì)模型，并基于分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)理論，推導(dǎo)了對(duì)應(yīng)的本構(gòu)關(guān)系和波動(dòng)方程.

本文基于Xu和Zhou(2014)提出的VTI介質(zhì)qP波方程，結(jié)合常Q模型的松弛函數(shù)，推導(dǎo)出新的黏聲各向異性純qP波方程，該方程沒有中間記憶變量，可以有效減少計(jì)算過程中的內(nèi)存存儲(chǔ)；該方程通過分?jǐn)?shù)階拉普拉斯算子將振幅衰減項(xiàng)和相位頻散項(xiàng)進(jìn)行解耦，在成像過程中可以有效進(jìn)行衰減補(bǔ)償；均勻TTI模型的數(shù)值模擬結(jié)果表明本文提出的黏聲各向異性波動(dòng)方程無偽橫波假象，可以穩(wěn)定高效地進(jìn)行數(shù)值求解.本文通過兩個(gè)VTI模型及一個(gè)TTI模型進(jìn)行逆時(shí)偏移成像測(cè)試，成像結(jié)果表明新的方程可以有效提高黏聲各向異性偏移成像結(jié)果的精度及分辨率.

1 黏聲各向異性純qP波方程

1.1 qP各向異性波動(dòng)方程

根據(jù)Tsvankin(1997)提出的VTI介質(zhì)中qP波的頻散關(guān)系：

(1)

其中，V(α)表示相對(duì)于對(duì)稱軸相位角為α的qP波相速度，VP0是qP波沿對(duì)稱軸的相速度，ε、δ均為Thomsen(1986)各向異性參數(shù)，f定義為

(2)

其中，VS0是qSV波沿對(duì)稱軸的相速度.根據(jù)qP波聲學(xué)近似假設(shè)VS0=0，且式(1)中所有參數(shù)都是隨空間變化的.

在式(1)及Alkhalifah(2000)的研究基礎(chǔ)上，解耦的qP波偏微分方程可以表示為

(3)

(4)

式(3)表示qP波在VTI介質(zhì)中的傳播.經(jīng)過變換，可以得到解耦的VTI介質(zhì)中的純qP波動(dòng)方程(Xu and Zhou，2014)：

(5)

其中，

(6)

考慮到式(5)中qP波是解耦的，我們將黏滯性引入到其中，以期得到純qP波在TI介質(zhì)中的黏聲波動(dòng)方程.

1.2 黏聲各向異性qP波動(dòng)方程

在廣義各向異性介質(zhì)中，常Q衰減模型的松弛函數(shù)可以表示為

(7)

其中，Γ(·)表示歐拉伽馬函數(shù)；Gmn表示介質(zhì)屬性的松弛函數(shù)；Cmn表示彈性剛度張量；βmn=arctan(1/Qmn)/π，是一個(gè)無量綱參數(shù)，且對(duì)于任意的品質(zhì)因子均有0<βmn<0.5；t0=1/ω0表示參考時(shí)間，ω0表示參考頻率.

基于歐拉伽馬函數(shù)的特性，各向異性黏彈介質(zhì)的廣義連續(xù)性方程可以表示為

σmn=Dmnemn,

(8)

(9)

對(duì)于TI介質(zhì)，矩陣Dmn具有與正交介質(zhì)相同的空元素，但只有5個(gè)獨(dú)立元素：D11，D13，D33，D55，D66.為了更加方便地描述VTI介質(zhì)中的地震傳播特征，我們可以使用Thomsen參數(shù)代替矩陣Dmn中的5個(gè)獨(dú)立元素，其關(guān)系式可以表示為

(10a)

(10b)

(10c)

(10d)

在這里，假定衰減是各向同性的，則參數(shù)βmn可以簡(jiǎn)化為β.因此，結(jié)合(8)—(10)三式，以及柯西方程和運(yùn)動(dòng)方程，Kjartansson模型中基于qP波頻散關(guān)系的黏聲VTI波動(dòng)方程可以表示為

(11)

當(dāng)β→0(Q→∞)時(shí)，上式可以退化為純qP波各向異性波動(dòng)方程；當(dāng)β→0.5(Q→0)時(shí)，上式可以退化為各向異性擴(kuò)散方程.為了避免存儲(chǔ)所有歷史時(shí)刻的波場(chǎng)，我們將分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)轉(zhuǎn)化為分?jǐn)?shù)階偽拉普拉斯算子(Zhu and Harris，2014)，因此，近似常Q模型的黏聲VTI波動(dòng)方程可以表示為

(12)

近似常Q模型的黏聲TTI波動(dòng)方程可以通過旋轉(zhuǎn)VTI介質(zhì)的坐標(biāo)系得到，可表示為

(13)

其中，

(14)

(15)

(16)

(17)

(18)

(19)

θ和φ分別是TTI介質(zhì)對(duì)稱軸的傾斜角度和傾斜方位角.

基于近似常Q模型的黏聲各向異性波動(dòng)方程包含兩項(xiàng)分?jǐn)?shù)階拉普拉斯算子，可以分別描述各向異性黏滯性介質(zhì)中地震波的能量衰減及相位頻散，并且將振幅衰減及相位頻散解耦表示，方程右端第一項(xiàng)描述相位頻散，第二項(xiàng)描述振幅衰減.其中，僅振幅衰減的黏聲TTI波動(dòng)方程可表示為

(20)

僅相位頻散的黏聲TTI波動(dòng)方程可表示為

(21)

為了說明黏聲各向異性純qP波動(dòng)方程描述的各向異性衰減介質(zhì)中地震波的傳播特征，此處使用二維均勻TTI模型進(jìn)行正演模擬.在均勻TTI模型中，橫縱向距離均為10 km，網(wǎng)格間距為10 m，設(shè)介質(zhì)速度為3500 m·s-1，密度均勻，各向異性參數(shù)ε為0.3，δ為0.25，θ為0.3，品質(zhì)因子Q為30，炮點(diǎn)位于模型中心，時(shí)間采樣間隔為0.5 ms，地震記錄總長(zhǎng)度為1.0 s.采用主頻為30 Hz的Ricker子波作為震源函數(shù).

圖1為均勻TTI模型在不同方程中的1.0 s波場(chǎng)快照對(duì)比圖.其中，(1)為TTI純qP波動(dòng)方程的模擬結(jié)果，(2)為僅振幅衰減的黏聲TTI純qP波動(dòng)方程的模擬結(jié)果，(3)為僅相位頻散的黏聲TTI純qP波動(dòng)方程的模擬結(jié)果，(4)為黏聲TTI純qP波動(dòng)方程的模擬結(jié)果.對(duì)比圖1中均勻衰減模型在不同方程下的波場(chǎng)快照可以發(fā)現(xiàn)，與純qP波各向異性波場(chǎng)快照相比，振幅衰減的黏聲TTI波場(chǎng)快照僅呈現(xiàn)出明顯的振幅衰減，波場(chǎng)能量減弱；相位頻散的黏聲TTI波場(chǎng)快照僅呈現(xiàn)出相位的滯后，振幅保持不變；而黏聲各向異性的波場(chǎng)快照中同時(shí)呈現(xiàn)出振幅衰減及相位畸變，可以表現(xiàn)出地震波在遇到衰減介質(zhì)之后的傳播特征.這說明黏聲各向異性波動(dòng)方程可以描述各向異性衰減介質(zhì)中地震波的運(yùn)動(dòng)學(xué)及動(dòng)力學(xué)特征，同時(shí)表明方程中的振幅衰減項(xiàng)及相位頻散項(xiàng)可以完全解耦，便于在逆時(shí)偏移過程中進(jìn)行衰減補(bǔ)償.

圖1 均勻TTI模型在不同方程中的波場(chǎng)快照(1.0 s)對(duì)比Fig.1 Four snapshot partsat 1.0s of the homogeneous TTI model with different wave equations

圖2為均勻模型在不同方程下的同一位置的單道地震記錄對(duì)比.與各向異性純qP波地震記錄相比，黏聲各向異性地震記錄的振幅明顯降低，能量減弱；相位出現(xiàn)明顯滯后，表現(xiàn)出相位延遲現(xiàn)象.

圖2 均勻TTI模型在不同方程中的同一單道記錄(5000 m處)對(duì)比Fig.2 Traces from Fig.1 of the attenuating homogeneous TTI model at horizontal 5000 m

圖3為基于近似常Q模型黏聲TTI波動(dòng)方程在不同品質(zhì)因子Q值下的波場(chǎng)快照切片對(duì)比，圖4為不同Q值下的單道記錄對(duì)比圖.從圖中可以看出，隨著Q值的減小，波場(chǎng)的振幅衰減及相位畸變程度逐漸增強(qiáng)，這說明文中推導(dǎo)的方程可以描述不同Q值下的地震波傳播特征.黏對(duì)于TTI介質(zhì)，傾斜角的引入使得直接采用偽譜法求解分?jǐn)?shù)階拉普拉斯算子時(shí)存在困難，因此在數(shù)值求解時(shí)需要先使用偽譜法求解空間偏導(dǎo)數(shù)，再將結(jié)果代入到兩個(gè)分?jǐn)?shù)階拉普拉斯算子中求取每個(gè)網(wǎng)格點(diǎn)處的值.

圖3 均勻TTI模型在不同Q值下的波場(chǎng)快照(1.0 s)Fig.3 Four snapshot parts at 1.0 s of the homogeneous TTI model with different Q values

1.3 黏聲各向異性純qP波逆時(shí)偏移

逆時(shí)偏移成像的理論基礎(chǔ)是Claerbout(1971)提出的時(shí)間一致性原理，其過程可以分為三個(gè)部分：震源波場(chǎng)的正向延拓，檢波點(diǎn)波場(chǎng)的逆向延拓，以及應(yīng)用成像條件進(jìn)行成像.在逆時(shí)偏移過程中，波場(chǎng)外推的目的是重建反射點(diǎn)處的入射和反射波場(chǎng).在黏滯性介質(zhì)中使用互相關(guān)成像條件進(jìn)行逆時(shí)偏移成像時(shí)，需要同時(shí)對(duì)震源波場(chǎng)和檢波點(diǎn)波場(chǎng)進(jìn)行衰減補(bǔ)償，使波場(chǎng)的能量得到恢復(fù)，相位畸變得到校正.在進(jìn)行波場(chǎng)外推時(shí)，Treeby和Cox(2010)及Zhu等(2014)指出，衰減補(bǔ)償可以通過反轉(zhuǎn)吸收衰減項(xiàng)的符號(hào)來實(shí)現(xiàn)，但保持等效頻散參數(shù)不變.因此，需要將振幅衰減項(xiàng)的符號(hào)由“+”號(hào)變換為“-”號(hào)來進(jìn)行振幅衰減的補(bǔ)償；同時(shí)，由于相位頻散項(xiàng)與時(shí)間延拓?zé)o關(guān)(即與頻率相關(guān)的相位速度在時(shí)間上保持不變)，保持相位頻散項(xiàng)的符號(hào)不變.

圖4 均勻TTI模型在不同Q值下的波場(chǎng)快照(1.0 s)單道對(duì)比Fig.4 Single traces from 1.0 s snapshots of homogeneous TTI model with different Q values

利用式(22)即可求解VTI介質(zhì)中Q補(bǔ)償?shù)恼鹪床▓?chǎng)Us：

(22)

同理，VTI介質(zhì)中Q補(bǔ)償?shù)臋z波點(diǎn)波場(chǎng)Ur可利用式(23)進(jìn)行求解：

(23)

在偏移成像的實(shí)際應(yīng)用中，地震資料中的較高頻率成分比較容易受到噪聲的影響.而在黏滯性介質(zhì)中，由于在波場(chǎng)逆時(shí)延拓中進(jìn)行了衰減補(bǔ)償，地震波能量呈指數(shù)增長(zhǎng)，可能會(huì)放大不需要的頻率成分，增加高頻噪聲，導(dǎo)致波場(chǎng)延拓的數(shù)值不穩(wěn)定.為了避免補(bǔ)償過程中高頻成分呈指數(shù)增長(zhǎng)而出現(xiàn)高波數(shù)發(fā)散問題，本文在計(jì)算補(bǔ)償波場(chǎng)時(shí)應(yīng)用低通濾波方法切除高頻噪聲，濾波器的截止波數(shù)根據(jù)介質(zhì)中最大速度上的截止頻率計(jì)算得到(Zhu et al.，2014).

2 數(shù)值算例

2.1 VTI-Hess模型逆時(shí)偏移

首先通過VTI-Hess模型進(jìn)行黏聲各向異性數(shù)值模擬，利用數(shù)值模擬得到的黏聲各向異性地震記錄分別進(jìn)行聲波、聲波各向異性、黏聲波及黏聲各向異性逆時(shí)偏移成像，對(duì)成像結(jié)果進(jìn)行對(duì)比，說明它們之間成像效果的差異.圖5為VTI-Hess模型，模型大小為9000 m×3750 m，空間間隔Δx=Δz=10 m.共布設(shè)了90炮，起始炮點(diǎn)位于(50 m,0 m)，炮間隔為100 m.每一炮均有601個(gè)接收點(diǎn)，平均分布在震源兩側(cè)，道間距為10 m.地震記錄時(shí)間長(zhǎng)度為6.0 s，時(shí)間采樣間隔為0.5 ms.采用主頻為30 Hz的Ricker子波作為震源函數(shù).

圖5 VTI-Hess模型(a) 速度模型； (b) ε模型； (c) δ模型； (d) Q模型.Fig.5 VTI-Hess model(a) Velocity model; (b) ε model; (c) δ model; (d) Q model.

為了驗(yàn)證本文所推導(dǎo)的黏聲各向異性波動(dòng)方程在描述地震波衰減特征方面的準(zhǔn)確性，分別進(jìn)行聲波各向異性及黏聲各向異性數(shù)值模擬，圖6為1.0 s及1.5 s時(shí)的波場(chǎng)快照對(duì)比.圖7為VTI-Hess模型的聲波各向異性單炮地震記錄及黏聲各向異性單炮地震記錄對(duì)比.圖8為使用不同方程模擬的單道(4500 m)地震記錄及對(duì)應(yīng)的振幅譜對(duì)比圖.

圖6 VTI-Hess模型數(shù)值模擬波場(chǎng)快照對(duì)比(a) 1.0 s純qP波VTI； (b) 1.0 s黏聲波VTI； (c) 1.5 s純qP波VTI； (d) 1.5 s黏聲波VTI.Fig.6 Snapshots of VTI-Hess model(a) Pure qP-wave VTI wave equation at 1.0 s; (b) Viscoacoustic VTI wave equation at 1.0 s; (c) Pure qP-wave VTI wave equation at 1.5 s; (d) Viscoacoustic VTI wave equation at 1.5 s.

圖7 VTI-Hess模型數(shù)值模擬地震記錄(a) 純qP波VTI波動(dòng)方程； (b) 黏聲VTI波動(dòng)方程.Fig.7 Records of VTI-Hess model(a) Record ofpure qP-wave VTI wave equation; (b) Record ofviscoacoustic VTI wave equation.

通過圖6及圖7中聲波與黏聲波數(shù)值模擬結(jié)果的對(duì)比可以發(fā)現(xiàn)，黏聲各向異性的數(shù)值模擬結(jié)果中呈現(xiàn)出顯著的振幅衰減及相位畸變，這說明文中推導(dǎo)的黏聲各向異性波動(dòng)方程在復(fù)雜介質(zhì)中可以較好地描述地震波的運(yùn)動(dòng)學(xué)及動(dòng)力學(xué)特征.從圖8中的單道(4500 m)地震記錄及對(duì)應(yīng)的振幅譜對(duì)比可以看出，在黏聲各向異性波動(dòng)方程中，由于地層的吸收衰減作用，地震資料的能量減弱，主頻降低，頻帶變窄，這會(huì)降低地震資料的分辨率.因此，必須對(duì)地震資料進(jìn)行吸收衰減補(bǔ)償逆時(shí)偏移.

圖8 同一單道(4500 m)地震記錄對(duì)比圖(a) 單道記錄對(duì)比圖； (b) 單道記錄振幅譜對(duì)比圖.Fig.8 Comparison of single traces with different equation(a) Comparison between single traces; (b) Comparison between amplitude spectra of single traces.

為了驗(yàn)證黏聲各向異性逆時(shí)偏移成像方法的效果，這里基于黏聲各向異性數(shù)值模擬結(jié)果進(jìn)行逆時(shí)偏移成像.圖9為VTI-Hess模型偏移結(jié)果.由圖可以看出，四種成像方法均可以對(duì)黏聲各向異性地震記錄進(jìn)行成像.進(jìn)一步對(duì)比可以發(fā)現(xiàn)，因未進(jìn)行衰減補(bǔ)償及忽略了各向異性，聲波各向同性逆時(shí)偏移結(jié)果(圖9a)中存在大量的構(gòu)造假象，反射波難以準(zhǔn)確歸位，同相軸聚焦性不好，衰減區(qū)域難以成像，總體成像質(zhì)量較差；進(jìn)行了衰減補(bǔ)償?shù)酿ぢ暩飨虍愋阅?/p>

圖9 VTI-Hess模型偏移剖面(a) 聲波逆時(shí)偏移剖面； (b) 黏聲波逆時(shí)偏移剖面； (c) 各向異性聲波逆時(shí)偏移剖面； (d) 各向異性黏聲波逆時(shí)偏移剖面.Fig.9 Migration images of VTI-Hess model(a) Acoustic RTM; (b) Acoustic VTI RTM; (c) Viscoacoustic RTM; (d) Viscoscoustic VTI RTM.

時(shí)偏移結(jié)果(圖9d)中的反射界面正確歸位，繞射波得到較好收斂，斷層、尖滅等復(fù)雜構(gòu)造也得到精細(xì)成像，衰減區(qū)域成像清晰，分辨率更高，能量更均衡，成像質(zhì)量較好.成像結(jié)果表明考慮介質(zhì)各向異性及黏滯性影響可有效改善成像質(zhì)量，提高成像剖面分辨率.

2.2 VTI-BP氣云模型逆時(shí)偏移

為了進(jìn)一步驗(yàn)證黏聲各向異性逆時(shí)偏移成像方法在復(fù)雜模型及強(qiáng)烈衰減區(qū)域的適用性，本節(jié)利用具有強(qiáng)烈衰減氣云構(gòu)造的VTI-BP氣云模型進(jìn)行數(shù)值模擬及逆時(shí)偏移試驗(yàn).圖10為VTI-BP氣云模型，模型大小為3970 m×1600 m，空間間隔Δx=Δz=10 m.共布設(shè)了79炮，起始炮點(diǎn)位于(35 m,10 m)，炮間隔為50 m.每一炮有401個(gè)接收點(diǎn)，平均分布在震源兩側(cè)，道間距為10 m.記錄時(shí)間長(zhǎng)度為2.0 s，時(shí)間采樣間隔為0.5 ms.采用主頻為30 Hz的Ricker子波作為震源函數(shù).

圖10 VTI-BP氣云模型(a) 速度模型； (b) ε模型； (c) δ模型； (d) Q模型.Fig.10 VTI-BP gas chimney model (a) Velocity model; (b) ε model; (c) δ model; (d) Q model.

圖11是基于VTI-BP氣云模型使用不同成像方法計(jì)算得到的0.5 s時(shí)震源波場(chǎng)及檢波點(diǎn)波場(chǎng)對(duì)比圖，其中左側(cè)一列為震源波場(chǎng)快照，右側(cè)一列為檢波點(diǎn)波場(chǎng)快照.圖11(a)和(e)中的波場(chǎng)是使用聲波逆時(shí)偏移算法得到的，均未進(jìn)行衰減補(bǔ)償，波場(chǎng)能量較弱，由于未考慮各向異性，波前面的位置不準(zhǔn)確；圖11(b)和(f)中的波場(chǎng)是使用黏聲各向同性逆時(shí)偏移算法得到的，經(jīng)過衰減補(bǔ)償?shù)牟▓?chǎng)能量得到恢復(fù)，但在未考慮各向異性的情況下波前面位置不準(zhǔn)確；圖11(c)和(g)中的波場(chǎng)是使用聲波各向異性逆時(shí)偏移算法得到的，圖中波前面的位置準(zhǔn)確，但因沒有考慮衰減補(bǔ)償?shù)挠绊憣?dǎo)致波場(chǎng)能量不均衡；圖11(d)和(h)中的波場(chǎng)是使用黏聲各向異性逆時(shí)偏移算法得到的，同時(shí)考慮了介質(zhì)各向異性及衰減補(bǔ)償因素的影響之后，波前面的位置準(zhǔn)確，波場(chǎng)能量得到恢復(fù).

圖11 VTI-BP氣云模型使用不同成像方法計(jì)算的震源波場(chǎng)及檢波點(diǎn)波場(chǎng)快照(0.5 s)對(duì)比(a) 聲波震源波場(chǎng)； (b) 聲波檢波點(diǎn)波場(chǎng)； (c) 黏聲震源波場(chǎng)； (d) 黏聲檢波點(diǎn)波場(chǎng)； (e) 聲波各向異性震源波場(chǎng)； (f) 聲波各向異性檢波點(diǎn)波場(chǎng)； (g) 黏聲各向異性震源波場(chǎng)； (h) 黏聲各向異性檢波點(diǎn)波場(chǎng).Fig.11 Snapshots at 0.5 s of source and receiver wavefields computed by different RTM of VTI-BP gas chimney model(a) Acoustic source wavefield snapshot; (b) Acoustic receiver wavefield snapshot; (c) Viscoacoustic source wavefield snapshot; (d) Viscoscoustic receiver wavefield snapshot; (e) Acoustic VTI source wavefield snapshot; (f) Acoustic VTI receiver wavefield snapshot; (g) Viscoacoustic VTI source wavefield snapshot; (h) Viscoscoustic VTI receiver wavefield snapshot.

圖12是VTI-BP氣云模型的偏移結(jié)果.圖13為VTI-BP氣云模型偏移結(jié)果在相同位置的局部放大對(duì)比圖.在聲波各向同性逆時(shí)偏移剖面中，因?yàn)槲纯紤]各向異性及吸收衰減的影響，導(dǎo)致偏移剖面中存在大量成像假象，繞射波無法收斂，衰減區(qū)域難以成像，整體能量不均衡，同相軸較難聚焦.在黏聲各向同性逆時(shí)偏移成像剖面中，雖然衰減區(qū)域的能量得到了一定補(bǔ)償，整體能量更加均衡；但由于忽略了各向異性參數(shù)的影響，繞射波難以收斂，反射界面難以準(zhǔn)確歸位，剖面中存在成像噪聲，部分同相軸連續(xù)性較差.在聲波各向異性逆時(shí)偏移剖面中，繞射波得到收斂，反射界面位置也比較準(zhǔn)確，但由于未進(jìn)行衰減補(bǔ)償，衰減區(qū)域的能量沒有得到較好的恢復(fù)，深部反射界面的成像不夠清晰，分辨率較低.而在黏聲各向異性逆時(shí)偏移剖面中，繞射波得到充分收斂，反射界面得到準(zhǔn)確歸位，衰減區(qū)域能量得到恢復(fù)，同相軸的連續(xù)性和聚焦性更好，細(xì)小構(gòu)造得到精細(xì)刻畫，剖面的能量更加均衡，分辨率得到提高，成像質(zhì)量更好.因此，黏聲各向異性逆時(shí)偏移成像方法可以對(duì)衰減區(qū)域的能量進(jìn)行補(bǔ)償，對(duì)深部層位清晰成像，使繞射波得到歸位，成像質(zhì)量良好，偏移剖面分辨率更高.

圖12 VTI-BP氣云模型偏移剖面(a) 聲波逆時(shí)偏移剖面； (b) 黏聲波逆時(shí)偏移剖面； (c) 聲波各向異性逆時(shí)偏移剖面； (d) 黏聲各向異性逆時(shí)偏移剖面.Fig.12 Migration images of VTI-BP gas chimney model(a) Acoustic RTM; (b) Acoustic VTI RTM; (c) Viscoacoustic RTM; (d) Viscoscoustic VTI RTM.

圖13 VTI-BP氣云模型偏移剖面局部放大圖(a) 聲波逆時(shí)偏移剖面； (b) 黏聲波逆時(shí)偏移剖面； (c) 各向異性聲波逆時(shí)偏移剖面； (d) 各向異性黏聲波逆時(shí)偏移剖面.Fig.13 Zoomed views of VTI-BP gas chimney model RTM images in Fig.10(a) Acoustic RTM; (b) Acoustic VTI RTM; (c) Viscoacoustic RTM; (d) Viscoacoustic VTI RTM.

2.3 TTI鹽丘模型逆時(shí)偏移

通過一個(gè)鹽丘模型進(jìn)行TTI介質(zhì)中黏聲波的逆時(shí)偏移成像，從而證明所提出的黏聲各向異性逆時(shí)偏移方法在復(fù)雜介質(zhì)中的適用性及偏移效果.圖14為TTI鹽丘模型，大小為8000 m×3500 m，空間間隔Δx=Δz=10 m.共布設(shè)了141炮，起始炮點(diǎn)位于(500 m,10 m)，炮間隔為50 m.每一炮有601個(gè)接收點(diǎn)，平均分布在震源兩側(cè)，道間距為10 m.記錄時(shí)間長(zhǎng)度為3.5 s，時(shí)間采樣間隔為0.5 ms.采用主頻為30 Hz的Ricker子波作為震源函數(shù).

圖14 TTI鹽丘模型(a) 速度模型； (b) ε模型； (c) δ模型； (d) θ模型； (e) Q模型.Fig.14 TTI-salt model(a) Velocity model; (b) ε model; (c) δ model; (d) θ model; (e) Q model.

圖15為TTI鹽丘模型的偏移結(jié)果.圖16為TTI鹽丘模型偏移結(jié)果(圖15)在相同位置的局部放大對(duì)比圖.從圖中可以看出，因?yàn)槲纯紤]各向異性參數(shù)及衰減補(bǔ)償?shù)挠绊?，聲波各向同性逆時(shí)偏移剖面中存在大量構(gòu)造假象，反射界面位置不準(zhǔn)確，剖面中能量不均衡，同相軸能量難以聚焦，剖面質(zhì)量較低.黏聲各向同性逆時(shí)偏移成像剖面中，在進(jìn)行吸收衰減補(bǔ)償之后，能量得到恢復(fù)，振幅較均衡；但由于忽略了各向異性參數(shù)的影響，繞射波難以收斂，反射界面不能準(zhǔn)確歸位，剖面中也存在一些成像噪聲，部分同相軸連續(xù)性較差，分辨率較低.在聲波各向異性逆時(shí)偏移剖面中，考慮了介質(zhì)各向異性對(duì)偏移結(jié)果的影響，繞射波得到收斂，反射界面位置準(zhǔn)確；但由于未進(jìn)行衰減補(bǔ)償，剖面的能量沒有得到恢復(fù)，深部反射界面成像不夠清晰.而在黏聲各向異性逆時(shí)偏移剖面中，同時(shí)考慮了各向異性及吸收衰減補(bǔ)償，繞射波得到充分收斂，反射界面得到準(zhǔn)確歸位，能量得到恢復(fù)，同相軸的連續(xù)性和聚焦性更好，分辨率得到提高，成像質(zhì)量更好.因此，與其他成像方法相比，黏聲各向異性逆時(shí)偏移可以對(duì)鹽丘下方高速衰減區(qū)域清晰成像，無成像噪聲，有效提高成像剖面質(zhì)量.

圖15 TTI鹽丘模型偏移剖面(a) 聲波逆時(shí)偏移剖面； (b) 黏聲波逆時(shí)偏移剖面； (c) 聲波各向異性逆時(shí)偏移剖面； (d) 黏聲各向異性逆時(shí)偏移剖面.Fig.15 Migration images of TTI-salt model(a) Acoustic RTM; (b) AcousticTTI RTM; (c) Viscoacoustic RTM; (d) Viscoacoustic TTI RTM.

3 分析與討論

3.1 成像方法的有效性和適用性分析

為了說明本文提出的基于黏聲各向異性波動(dòng)方程成像方法的有效性和適用性，本節(jié)基于VTI-BP氣云模型，采用如下方程(Xu et al.，2015)進(jìn)行數(shù)值模擬：

(24)

在進(jìn)行數(shù)值模擬時(shí)，參數(shù)與2.2節(jié)中的參數(shù)保持一致.基于方程(24)的數(shù)值模擬結(jié)果，本節(jié)利用1.3節(jié)中的成像方法進(jìn)行衰減補(bǔ)償?shù)钠瞥上?，并與2.2節(jié)中的黏聲各向異性逆時(shí)偏移結(jié)果(圖12d)進(jìn)行對(duì)比.

圖16 TTI鹽丘模型偏移剖面局部放大圖(a) 聲波逆時(shí)偏移剖面； (b) 黏聲波逆時(shí)偏移剖面； (c) 聲波各向異性逆時(shí)偏移剖面； (d) 黏聲各向異性逆時(shí)偏移剖面.Fig.16 Zoomed views of TTI-salt model RTM images in Fig.13(a) Acoustic RTM; (b) Acoustic TTI RTM; (c) Viscoacoustic RTM; (d) Viscoacoustic TTI RTM.

對(duì)比圖17中的偏移結(jié)果，可以看到基于方程(24)模擬數(shù)據(jù)的偏移剖面與基于本文推導(dǎo)的黏聲各向異性波動(dòng)方程模擬數(shù)據(jù)的偏移剖面一致，深層衰減區(qū)域的能量得到了恢復(fù)，層位構(gòu)造清晰，證明了本文成像方法的有效性和適用性.

圖17 VTI-BP氣云模型黏聲各向異性逆時(shí)偏移剖面(a) 基于方程(12)模擬數(shù)據(jù)； (b) 基于方程(24)模擬數(shù)據(jù).Fig.17 Viscoacoustic VTI RTM images of VTI-BP gas chimney model using the synthetic data obtained with (a) equation (12) and (b) equation (24)

3.2 低通濾波器對(duì)偏移剖面質(zhì)量的影響

在黏性介質(zhì)中，由于需要在波場(chǎng)逆時(shí)延拓過程中對(duì)地震波吸收衰減進(jìn)行補(bǔ)償，地震波能量呈指數(shù)型增長(zhǎng)，可能會(huì)放大不需要的頻率成分，增加高頻噪聲，導(dǎo)致波場(chǎng)延拓的數(shù)值不穩(wěn)定.為了避免這個(gè)問題，本文應(yīng)用低通濾波方法切除高頻噪聲，濾波器的截止波數(shù)根據(jù)介質(zhì)中最大速度上的截止頻率計(jì)算得到.但這會(huì)使最終的偏移剖面損失一些有效的高波數(shù)成分，一定程度上會(huì)降低偏移剖面的精度.本節(jié)將對(duì)比不同濾波器參數(shù)下的偏移結(jié)果，以更直觀地說明低通濾波器對(duì)最終偏移剖面質(zhì)量的影響.

在現(xiàn)階段的衰減補(bǔ)償成像方法中，主要通過低通濾波器來保證波場(chǎng)延拓的穩(wěn)定性，濾波參數(shù)的選取大多是依據(jù)經(jīng)驗(yàn)進(jìn)行取值.從圖18中可以看出，當(dāng)濾波參數(shù)過小或過大時(shí)，都會(huì)導(dǎo)致偏移剖面質(zhì)量的降低.當(dāng)?shù)屯V波器的截止波數(shù)過小時(shí)，說明截去的高波數(shù)成分較多，會(huì)導(dǎo)致深層構(gòu)造成像不清晰，降低成像剖面質(zhì)量；當(dāng)?shù)屯V波器的截止波數(shù)過大時(shí)，表明截去的高波數(shù)成分較少，但無法保證波場(chǎng)延拓的穩(wěn)定性，降低偏移剖面的精度.因此，在使用低通濾波器時(shí)，截止波數(shù)的選取一定程度上會(huì)影響衰減補(bǔ)償偏移剖面的質(zhì)量，這要求我們?cè)诮窈蟮难芯恐袑ふ乙环N科學(xué)、定量的保持波場(chǎng)延拓穩(wěn)定的方法.

圖18 不同截止波數(shù)偏移剖面(a) 截止波數(shù)偏小(0.05)； (b) 截止波數(shù)適中(0.2)； (c) 截止波數(shù)偏大(0.35).Fig.18 Migration images of different cut-off wavenumbers(a) Low cut-off wavenumber (0.05); (b) Moderate cut-off wavenumber (0.2); (c) High cut-off wavenumber (0.35).

4 結(jié)論

(1) 本文在各向異性純qP波動(dòng)方程的基礎(chǔ)上，結(jié)合常Q模型的松弛函數(shù)，推導(dǎo)出黏聲各向異性純qP波動(dòng)方程.均勻模型數(shù)值模擬結(jié)果表明，該方程可以正確描述TI黏聲介質(zhì)中地震波的運(yùn)動(dòng)學(xué)及動(dòng)力學(xué)特征，并且方程中的振幅衰減項(xiàng)與相位頻散項(xiàng)解耦表示，便于在偏移過程中進(jìn)行衰減補(bǔ)償.

(2) 通過黏聲各向異性逆時(shí)偏移成像結(jié)果可以看出，同時(shí)考慮介質(zhì)各向異性及衰減補(bǔ)償?shù)某上窠Y(jié)果可以較好恢復(fù)地下衰減區(qū)域的能量，使反射界面準(zhǔn)確歸位，繞射波得到收斂，增強(qiáng)同相軸的連續(xù)性和聚焦性，精細(xì)刻畫細(xì)小的地質(zhì)構(gòu)造，剖面中能量得到均衡，提高剖面分辨率，成像質(zhì)量更好.

(3) 本文通過低通濾波來保證波場(chǎng)外推時(shí)的穩(wěn)定性，這會(huì)使最終的偏移剖面損失一些有效的高波數(shù)成分，一定程度上會(huì)降低偏移剖面的精度.因此，在今后的研究中可以進(jìn)一步探討上述問題.