結合泊松算法和波場分解成像條件的各向異性優化純聲波方程逆時偏移

徐世剛 包乾宗 任志明 劉 洋

(①長安大學地質工程與測繪學院地球物理系，陜西西安 710054； ②中國石油大學(北京)油氣資源與探測國家重點實驗室，北京 102249)

0 引言

地下介質普遍呈現各向異性，橫向各向同性(TI)介質是地震勘探學界研究最廣的一類各向異性模型，根據其對稱軸方向可以分為VTI、HTI和TTI模型。逆時偏移技術受傾角限制小，能夠適應復雜速度模型，獲得高精度成像剖面，因此各向異性逆時偏移研究具有重要意義。

相比于各向異性彈性波方程，各向異性聲波方程形式簡單、計算量較小、波場特征明確。目前在各向異性波場模擬和偏移中主要采用偽聲波方程，該類方程的一種推導思路是將耦合的各向異性P-SV波頻散關系中的橫波速度設置為0簡化而來[1]。在此基礎上，發展了多種形式的偽聲波方程[2-8]。另一種偽聲波方程的推導思路是基于聲學近似的胡克定律[9-10]。但是當給定的各向異性參數不滿足偽聲波假設時，現有偽聲波方程在進行P波模擬時，容易出現偽橫波干擾及數值不穩定。為了徹底解決偽聲波存在的問題，Liu等[11]解耦了VTI介質P-SV波頻散關系，獲得單獨的P波和SV波方程。相比于顯式表達的偽聲波方程，各向異性純聲波方程包含復雜的擬微分算子，常規差分法難以直接求解。迄今為止，已經發展了多種數值算法求解純聲波方程，如偽譜法/偽解析法[12-13]、低秩分解法[14]、迭代求解法[15]、混合求解法[16]和快速泊松算法[17-18]等。目前關于純聲波方程的研究主要集中于平衡計算精度和效率方面，以及在偏移成像中的進一步應用等[19-23]。

傳統逆時偏移技術對震源波場和檢波點波場做互相關，結果中容易出現較強的低頻噪聲。Zhang等[24]采用拉普拉斯濾波壓制低頻噪聲； Yoon等[25]、康瑋等[26]利用Poynting矢量切除大角度成像結果，從而壓制低頻噪聲； Liu等[27]提出波場分解構建逆時偏移成像條件，選擇不同波場做互相關，能夠較好地抑制成像噪聲。在Liu等[27]的研究基礎上，基于Hilbert變換構建復數波場的策略被用于提高計算效率[28-32]。波場分解成像條件已應用于各向同性聲波和彈性波逆時偏移，考慮到該方法的優勢，本文將其推廣到各向異性介質逆時偏移。

作為對比，本文首先回顧了兩種常用的TTI偽聲波方程，分析了其存在的問題； 然后利用最小二乘方法求取了TTI優化純聲波頻散關系，并采用泊松算法和有限差分相結合的數值方案進行求解； 最后結合介質各向異性特征對波場分解成像條件進行改進，將其推廣至純聲波逆時偏移。相關理論和算例表明，將優化純聲波方程和波場分解成像條件聯合應用于各向異性純聲波逆時偏移，能有效壓制偽橫波及低頻噪聲，獲得高質量成像結果。

1 理論方法

1.1 基于限制性橫波速度的TTI偽聲波方程

通過將原始P-SV波頻散關系中的橫波速度設置為0，Alkhalifah[1]推導了各向異性偽聲波方程，該方程能夠有效簡化計算，其多種改進版本被廣泛應用于各向異性波場模擬和偏移成像[2-8]。但直接將垂向橫波速度設置為0并不能完全滿足各向異性介質中的物理學規律，當各向異性參數不滿足近似條件時(即橢圓近似條件)，偽聲波方程模擬容易出現SV波干擾及不穩定。為此，學者們提出了不同的解決方案，其中一種常用的基于限制性橫波速度的TTI偽聲波方程形式[4]為

(1)

(2)

其中θ為對稱軸傾角。特別地，通過引入限制性橫波速度能夠有效提高穩定性。當vSz=0時，式(1)退化為傳統偽聲波方程，形式為

(3)

采用有限差分即可求解式(1)和式(3)。圖1為兩個方程計算的雙層TTI模型(表1)的波場切片，其中上、下層介質厚度均為1500m。

表1 雙層TTI模型參數

由圖1可以看出：當介質具有橢圓各向異性特征(ε=δ)時，兩種方程均可產生高精度P波響應，如圖1左所示； 當各向異性參數不滿足橢圓假設時(ε≠δ)，兩種方程會產生較強的偽橫波污染，如圖1中和圖1a右所示； 當ε<δ時，橫波速度為0的偽聲波方程(式(3))模擬出現不穩定，如圖1b右所示。

1.2 TTI優化純聲波方程及其結合泊松算法和有限 差分的解法

多種改進版本的各向異性偽聲波方程僅能在一定程度上壓制偽橫波干擾，提高穩定性，但不能從根本上解決問題。眾多研究表明，發展純聲波方程能夠徹底消除偽橫波污染及不穩定。相比于顯式表達的偽聲波方程，純聲波方程包含復雜的偏微分算子，常規差分法難以直接求解。考慮到計算精度和效率，本文給出一種改進的、結合泊松算法和有限差分的各向異性優化純聲波方程求解方法[17-18]。

解耦的TTI純聲波頻散關系[11]為

(4)

(5)

將式(5)代入式(4)，可得到目前最常用的各向異性純聲波頻散關系。考慮到一階泰勒展開的精度較低，高階泰勒展開的計算量較大，為了平衡純聲波模擬時的精度和效率，本文借鑒優化的思想[17-18]獲得一種高精度的純聲波頻散關系。對式(4)中的根號項做如下近似[17-18]

圖1 雙層TTI模型式(1)(a)和式(3)方程模擬的波場切片(b)左：模型1； 中：模型2； 右：模型3

(6)

式(6)相對于待求近似系數(f1，f2，f3)為線性方程，因此可以通過極小化兩端的誤差獲得近似系數值。構建的目標函數為

(7)

通過最小二乘方法求解上式，即可獲得優化系數(f1，f2，f3)。

將式(6)代入式(4)，整理可得優化的純聲波頻散關系

(8)

式中

(9)

特別地，當w1=1+2ε、w2=1、w3=-2(ε-δ)時，式(8)退化為常用的一階泰勒展開純聲波頻散式。式(8)在時間—空間域可寫為[17]

(10)

引入輔助波場Q，將式(10)分解為[15，17]

(11)

(12)

式(11)為顯式方程，采用常規差分即可求解。式(12)為隱式方程，其計算過程涉及線性方程組求解。Chu等[15]采用線性迭代求解隱式方程，但計算量較大。為了提高效率，本文給出一種基于泊松算法的實現策略[17]。

將式(11)和式(12)進一步整理為

(13)

?2Q=P

(14)

式(13)和式(14)組成了新的TTI優化純聲波方程，其中式(13)可采用常規差分法求解。式(14)為泊松方程，基于前人研究成果，可通過泊松算法實現，計算步驟[17]如下。

(1)采用五點差分離散式(14)

Qi-1,j+Qi+1,j+Qi,j-1+Qi,j+1-4Qi,j=h2Pi,j

(15)

式中：h為網格間距；i、j分別為x、z方向網格點序號。

(16)

式中：Lx和Lz分別為x、z方向網格點數；n、l分別為波數域兩個方向的網格點序號。

(3)計算輔助波場Q的傅里葉響應

(17)

(4)求解輔助波場Q

(18)

通過顯式差分離散式(13)，通過式(15)～式(18)可實現式(14)的求解。這樣能夠保證純聲波方程式(13)和式(14)的高效求解。

圖2對比了一階泰勒展開方程(式(5))和優化展開方程(式(6))對根號部分的近似精度。由圖可見，與精確值(圖2a)相比，一階泰勒展開存在較大誤差(圖2b)，而最小二乘優化展開能夠獲得更高的近似精度(圖2c)。圖3為優化純聲波方程模擬的雙層TTI模型波場切片，與圖1偽聲波方程的計算結果相比，純聲波方程在不同參數條件下均能產生精確有效的P波響應，同時能夠完全消除偽橫波和數值不穩定。

1.3 改進波場分解成像條件

在逆時偏移成像中，互相關成像條件具有實現簡單、計算高效等優點，得到廣泛應用。常規互相關成像條件[34]為

(19)

式中：x為空間坐標；I(x)為成像結果；Ps(x,t)為震源波場；Pr(x,t)為檢波點波場；tmax為記錄長度。

常規互相關成像條件易引起較強低頻噪聲，原因是互相關過程中沿所有方向傳播的震源波場和檢波點波場均參與成像，同向傳播的波場之間互相關疊加，易形成較強的低頻噪聲。為了壓制傳統互相關成像過程中的低頻噪聲，除了對低頻噪聲進行濾波處理外[24]，另一種有效策略是對參與成像的全波場沿不同方向進行分解，選取傳播方向不同的分量波場參與成像，能夠有效減弱低頻噪聲干擾[27-32]。

以二維情況為例，對式(19)中的震源波場Ps(x,t)和檢波點波場Pr(x,t)進行分解

(20)

(21)

式中上標“↑”、“↓”、“←”、“→”分別表示上行、下行、左行和右行波場分量。

由式(19)～式(21)可知，震源波場分量與檢波點波場分量之間兩兩互相關可得16種成像剖面。其中傳播方向相反的波場分量互相關能夠產生有效的成像結果，而同向傳播的波場分量之間互相關主要產生低頻噪聲[27-32]。在前人關于波場分解成像條件的研究基礎上，本文采用一種基于Hilbert變換構建復數波場的分解算法。該分解方法能夠進行全波場的顯式分離，目前被廣泛用于各向同性介質的

圖2 均勻各向異性模型的精確微分算子(a)及一階泰勒展開誤差(b)與最小二乘優化展開誤差(c)的對比模型參數為ε=0.2，δ=0.1，θ=60°

圖3 雙層TTI模型優化純聲波方程模擬的波場切片(a)模型1； (b)模型2； (c)模型3

逆時偏移[27-32]。本文將其推廣到各向異性介質的逆時偏移，實現步驟如下。

(1)基于式(13)和式(14)，震源波場的復數外推形式可以表示為

(22)

(2)基于步驟(1)構建的復數波場，以震源波場為例，上行波場可以表示為

(23)

(24)

圖4 雙層TTI模型4的波場切片(a)及分離的上(b)、下(c)、左(d)、右行波場(e)

2 模型算例

2.1 地塹模型

各向異性地塹模型的尺寸為4000m×3000m，相關參數如圖5所示。網格間距為10m，時間步長為1ms，震源選用主頻為28Hz的Ricker子波。

考慮到橫波速度為0的偽聲波方程存在穩定性問題，圖6為基于限制性橫波速度的偽聲波和優化純聲波方程模擬的波場切片，可以看出，相比于偽聲波方程的模擬結果(圖6a)，純聲波方程能夠有效壓

圖5 各向異性地塹模型

制偽橫波干擾，產生高精度P波響應(圖6b)。

圖7為圖6b的純聲波波場沿不同方向分解結果，可以發現，上、下、左、右行波得到了有效分離。傳統互相關成像條件和波場分解成像條件的偏移結果如圖8所示。由圖8a可知，基于傳統互相關條件的成像結果容易產生較強低頻噪聲干擾。對比部分分量波場的成像結果可見同向傳播的波場之間互相關主要產生低頻噪聲，波場分解成像條件能夠較好地將這部分干擾剔除，反向傳播波場之間互相關能夠產生有效的成像剖面。對圖8a中的成像剖面進行濾波處理，結果如圖8b所示。由圖可知，不同剖面的低頻噪音得到較好地去除。反向傳播波場成像

圖6 地塹模型不同方程模擬的波場切片(a)偽聲波方程(式(1))； (b)優化純聲波方程(式(13)、式(14))

圖7 地塹模型分離的上(a)、下(b)、左(c)、右行波場(d)

圖8 地塹模型優化純聲波方程不同條件成像結果(a)及其拉普拉斯濾波結果(b)左：全波場互相關； 中：同向傳播波場互相關之和； 右：反向傳播波場互相關之和

結果之和略優于傳統全波場成像結果。相關結果表明，不同方向行波之間的互相關結果可以作為傳統方法的補充，成像結果精度更高。

2.2 修改的Marmousi TTI模型

修改的Marmousi TTI模型，如圖9所示。將計算區域剖分為500×350個網格單元。用于逆時偏移的速度場和各向異性參數場均基于五點平滑獲得。圖10為各向同性聲波、各向異性偽聲波和優化純聲波的成像結果，相比之下，基于優化純聲波方程的逆時偏移能夠獲得更高精度的成像。圖11為基于優化純聲波方程的不同條件成像結果的對比，可見，傳統全波場參與的互相關成像，淺層存在嚴重的低頻噪聲污染。同向傳播的波場之間互相關主要產生低頻成像噪聲，波場分解成像條件能夠較好地將低頻噪聲去除，反向傳播波場之間互相關能夠產生高精度的成像結果。相比于經過拉普拉斯濾波的成像剖面而言，基于波場分解成像條件的偏移結果不需要濾波處理，能夠較好地保證振幅信息。總體而言，基于波場分解的成像方法可以作為傳統方法的有效補充，以獲得不同展布方向地下構造的高精度成像。

圖9 修改的Marmousi TTI模型(a)vPz； (b)ε； (c)δ； (d)θ

圖10 修改的Marmousi TTI模型不同方程的成像結果(拉普拉斯濾波后)(a)各向同性聲波方程； (b)偽聲波方程(式(1))；(c)優化純聲波方程(式(13)、式(14))

圖11 修改的Marmousi TTI模型的優化純聲波方程不同條件成像結果

3 結束語

本文首先回顧了TTI介質中兩種常用的偽聲波方程，在各向異性參數不滿足聲學近似時，利用該方程模擬容易產生偽橫波干擾及不穩定現象。為了解決現有偽聲波方程存在的問題，應用最小二乘方法獲得了優化的純聲波頻散關系，在此基礎上采用一種泊松算法和有限差分相結合的高效精確各向異性純聲波求解算法，能夠有效避免偽橫波干擾及數值不穩定。將各向同性介質中基于Hilbert變換的復數波場分解成像條件推廣到TTI介質優化純聲波逆時偏移成像，選擇傳播方向相反的分量波場參與最終成像。理論和數值算例驗證了本文方法能夠有效壓制各向異性逆時偏移的偽橫波干擾和低頻噪聲，獲得高精度的成像結果。