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

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

中國石油大學(華東)深層油氣重點實驗室， 青島 266580

0 引言

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

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

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

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

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

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

1.1 qP各向異性波動方程

根據Tsvankin(1997)提出的VTI介質中qP波的頻散關系：

(1)

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

(2)

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

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

(3)

(4)

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

(5)

其中，

(6)

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

1.2 黏聲各向異性qP波動方程

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

(7)

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

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

σmn=Dmnemn,

(8)

(9)

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

(10a)

(10b)

(10c)

(10d)

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

(11)

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

(12)

近似常Q模型的黏聲TTI波動方程可以通過旋轉VTI介質的坐標系得到，可表示為

(13)

其中，

(14)

(15)

(16)

(17)

(18)

(19)

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

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

(20)

僅相位頻散的黏聲TTI波動方程可表示為

(21)

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

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

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

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

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

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

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

1.3 黏聲各向異性純qP波逆時偏移

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

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

利用式(22)即可求解VTI介質中Q補償的震源波場Us：

(22)

同理，VTI介質中Q補償的檢波點波場Ur可利用式(23)進行求解：

(23)

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

2 數值算例

2.1 VTI-Hess模型逆時偏移

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

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

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

圖6 VTI-Hess模型數值模擬波場快照對比(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模型數值模擬地震記錄(a) 純qP波VTI波動方程； (b) 黏聲VTI波動方程.Fig.7 Records of VTI-Hess model(a) Record ofpure qP-wave VTI wave equation; (b) Record ofviscoacoustic VTI wave equation.

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

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

為了驗證黏聲各向異性逆時偏移成像方法的效果，這里基于黏聲各向異性數值模擬結果進行逆時偏移成像.圖9為VTI-Hess模型偏移結果.由圖可以看出，四種成像方法均可以對黏聲各向異性地震記錄進行成像.進一步對比可以發現，因未進行衰減補償及忽略了各向異性，聲波各向同性逆時偏移結果(圖9a)中存在大量的構造假象，反射波難以準確歸位，同相軸聚焦性不好，衰減區域難以成像，總體成像質量較差；進行了衰減補償的黏聲各向異性逆

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

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

2.2 VTI-BP氣云模型逆時偏移

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

圖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氣云模型使用不同成像方法計算得到的0.5 s時震源波場及檢波點波場對比圖，其中左側一列為震源波場快照，右側一列為檢波點波場快照.圖11(a)和(e)中的波場是使用聲波逆時偏移算法得到的，均未進行衰減補償，波場能量較弱，由于未考慮各向異性，波前面的位置不準確；圖11(b)和(f)中的波場是使用黏聲各向同性逆時偏移算法得到的，經過衰減補償的波場能量得到恢復，但在未考慮各向異性的情況下波前面位置不準確；圖11(c)和(g)中的波場是使用聲波各向異性逆時偏移算法得到的，圖中波前面的位置準確，但因沒有考慮衰減補償的影響導致波場能量不均衡；圖11(d)和(h)中的波場是使用黏聲各向異性逆時偏移算法得到的，同時考慮了介質各向異性及衰減補償因素的影響之后，波前面的位置準確，波場能量得到恢復.

圖11 VTI-BP氣云模型使用不同成像方法計算的震源波場及檢波點波場快照(0.5 s)對比(a) 聲波震源波場； (b) 聲波檢波點波場； (c) 黏聲震源波場； (d) 黏聲檢波點波場； (e) 聲波各向異性震源波場； (f) 聲波各向異性檢波點波場； (g) 黏聲各向異性震源波場； (h) 黏聲各向異性檢波點波場.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氣云模型的偏移結果.圖13為VTI-BP氣云模型偏移結果在相同位置的局部放大對比圖.在聲波各向同性逆時偏移剖面中，因為未考慮各向異性及吸收衰減的影響，導致偏移剖面中存在大量成像假象，繞射波無法收斂，衰減區域難以成像，整體能量不均衡，同相軸較難聚焦.在黏聲各向同性逆時偏移成像剖面中，雖然衰減區域的能量得到了一定補償，整體能量更加均衡；但由于忽略了各向異性參數的影響，繞射波難以收斂，反射界面難以準確歸位，剖面中存在成像噪聲，部分同相軸連續性較差.在聲波各向異性逆時偏移剖面中，繞射波得到收斂，反射界面位置也比較準確，但由于未進行衰減補償，衰減區域的能量沒有得到較好的恢復，深部反射界面的成像不夠清晰，分辨率較低.而在黏聲各向異性逆時偏移剖面中，繞射波得到充分收斂，反射界面得到準確歸位，衰減區域能量得到恢復，同相軸的連續性和聚焦性更好，細小構造得到精細刻畫，剖面的能量更加均衡，分辨率得到提高，成像質量更好.因此，黏聲各向異性逆時偏移成像方法可以對衰減區域的能量進行補償，對深部層位清晰成像，使繞射波得到歸位，成像質量良好，偏移剖面分辨率更高.

圖12 VTI-BP氣云模型偏移剖面(a) 聲波逆時偏移剖面； (b) 黏聲波逆時偏移剖面； (c) 聲波各向異性逆時偏移剖面； (d) 黏聲各向異性逆時偏移剖面.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) 聲波逆時偏移剖面； (b) 黏聲波逆時偏移剖面； (c) 各向異性聲波逆時偏移剖面； (d) 各向異性黏聲波逆時偏移剖面.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鹽丘模型逆時偏移

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

圖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鹽丘模型的偏移結果.圖16為TTI鹽丘模型偏移結果(圖15)在相同位置的局部放大對比圖.從圖中可以看出，因為未考慮各向異性參數及衰減補償的影響，聲波各向同性逆時偏移剖面中存在大量構造假象，反射界面位置不準確，剖面中能量不均衡，同相軸能量難以聚焦，剖面質量較低.黏聲各向同性逆時偏移成像剖面中，在進行吸收衰減補償之后，能量得到恢復，振幅較均衡；但由于忽略了各向異性參數的影響，繞射波難以收斂，反射界面不能準確歸位，剖面中也存在一些成像噪聲，部分同相軸連續性較差，分辨率較低.在聲波各向異性逆時偏移剖面中，考慮了介質各向異性對偏移結果的影響，繞射波得到收斂，反射界面位置準確；但由于未進行衰減補償，剖面的能量沒有得到恢復，深部反射界面成像不夠清晰.而在黏聲各向異性逆時偏移剖面中，同時考慮了各向異性及吸收衰減補償，繞射波得到充分收斂，反射界面得到準確歸位，能量得到恢復，同相軸的連續性和聚焦性更好，分辨率得到提高，成像質量更好.因此，與其他成像方法相比，黏聲各向異性逆時偏移可以對鹽丘下方高速衰減區域清晰成像，無成像噪聲，有效提高成像剖面質量.

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

3 分析與討論

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

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

(24)

在進行數值模擬時，參數與2.2節中的參數保持一致.基于方程(24)的數值模擬結果，本節利用1.3節中的成像方法進行衰減補償的偏移成像，并與2.2節中的黏聲各向異性逆時偏移結果(圖12d)進行對比.

圖16 TTI鹽丘模型偏移剖面局部放大圖(a) 聲波逆時偏移剖面； (b) 黏聲波逆時偏移剖面； (c) 聲波各向異性逆時偏移剖面； (d) 黏聲各向異性逆時偏移剖面.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.

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

圖17 VTI-BP氣云模型黏聲各向異性逆時偏移剖面(a) 基于方程(12)模擬數據； (b) 基于方程(24)模擬數據.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 低通濾波器對偏移剖面質量的影響

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

在現階段的衰減補償成像方法中，主要通過低通濾波器來保證波場延拓的穩定性，濾波參數的選取大多是依據經驗進行取值.從圖18中可以看出，當濾波參數過小或過大時，都會導致偏移剖面質量的降低.當低通濾波器的截止波數過小時，說明截去的高波數成分較多，會導致深層構造成像不清晰，降低成像剖面質量；當低通濾波器的截止波數過大時，表明截去的高波數成分較少，但無法保證波場延拓的穩定性，降低偏移剖面的精度.因此，在使用低通濾波器時，截止波數的選取一定程度上會影響衰減補償偏移剖面的質量，這要求我們在今后的研究中尋找一種科學、定量的保持波場延拓穩定的方法.

圖18 不同截止波數偏移剖面(a) 截止波數偏小(0.05)； (b) 截止波數適中(0.2)； (c) 截止波數偏大(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 結論

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

(2) 通過黏聲各向異性逆時偏移成像結果可以看出，同時考慮介質各向異性及衰減補償的成像結果可以較好恢復地下衰減區域的能量，使反射界面準確歸位，繞射波得到收斂，增強同相軸的連續性和聚焦性，精細刻畫細小的地質構造，剖面中能量得到均衡，提高剖面分辨率，成像質量更好.

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